	program rf6uP
c
c	Copyright 2002-2003  Leo Breiman and Adele Cutler
c
c	This is free open source software but its use, in part or 
c	in whole, in any commercial product that is sold for profit 
c	is prohibited without the written consent of Leo Breiman 
c	and Adele Cutler.
c
c	We very much appreciate bug notices and suggested improvements.
c
c	leo@stat.berkeley.edu   adele@math.usu.edu
c
c
c	SET ALL PARAMETERS FIRST GROUP BELOW.  GENERALLY, 
c	SETTING PARAMETERS TO ZERO TURNS THE CORRESPONDING
c	OPTION OFF. 
c
c	ALL RELEVANT OUTPUT FILES MUST BE GIVEN NAMES--SEE BELOW.
c	NOTE to me--turning on classwt produces a lot of empty nodes.
c
c	Note to adele: get rid of the following declarations for distribution!
	integer mdim, nsample0, nclass, maxcat, ntest, labelts, labeltr,
     &	mtry0, ndsize, jbt, look, lookcls, jclasswt, mdim2nd, mselect,
     &	imp, interact, impn, impfast, nprdata, nprox, nrnn, noutlier,
     &	nscale, nprot, mfixall, mfixtr, nsample, nrnodes,
     &	mimp, nimp, near, nearts, ntest0, nscale0, nprot0, 
     &	isaverf, isavepar, isavefill, isaveprox, 
     &	irunrf, ireadpar, ireadfill, ireadprox,
     &	isumout, idataout, impfastout, impout, impnout, interout, 
     &	iprotout, iproxout, iscaleout, ioutlierout 
	real code
c
c	-------------------------------------------------------
c	CONTROL PARAMETERS
c
	parameter(
c  		DESCRIBE DATA		
     1		mdim=36, nsample0=4435, nclass=6, maxcat=1,
     1		ntest=0, labelts=0, labeltr=1,
c
c		SET RUN PARAMETERS
     2		mtry0=6, ndsize=1, jbt=50, look=10, lookcls=1,
c
c		SET RUN OPTIONS
     3		jclasswt=0, mdim2nd=0, mselect=0,
c
c		SET IMPORTANCE OPTIONS	
     4		imp=0, interact=0, impn=0, impfast=0,
c
c		SET PROXIMITY COMPUTATIONS
     5		nprox=0, nrnn=5,
c
c		SET OPTIONS BASED ON PROXIMITIES
     6		noutlier=0, nscale=0, nprot=0,
c
c		REPLACE MISSING VALUES  
     7		code=-999, mfixall=0, mfixtr=0, 
c
c		SAVING A FOREST
     8		isaverf=1, isavepar=0, isavefill=0, isaveprox=0,
c
c		RUNNING A SAVED FOREST
     9		irunrf=0, ireadpar=0, ireadfill=0, ireadprox=0)
c
c	-------------------------------------------------------
c	OUTPUT CONTROLS	
c
	parameter(
     &		isumout	=  	1, 	!0/1	1=summary to screen
     &		idataout=	1,   	!0/1/2	1=train,2=adds test	(7)
     &		impfastout=	1,	!0/1	1=gini fastimp		(8)
     &		impout=		1,   	!0/1/2 	1=imp,2=to screen	(9)
     &		impnout=	1,   	!0/1	1=impn			(10)
     &		interout=	1,	!0/1/2	1=interaction,2=screen	(11)
     &		iprotout=	1,	!0/1/2	1=prototypes,2=screen	(12)
     &		iproxout=	1,	!0/1/2	1=prox,2=adds test	(13)
     &		iscaleout=	1,	!0/1	1=scaling coors		(14)
     &		ioutlierout=	1)	!0/1/2	1=train,2=adds test 	(15)
c
c	-------------------------------------------------------
c	DERIVED PARAMETERS (DO NOT CHANGE)
c
 	parameter(
     &	nsample=(2-labeltr)*nsample0,
     &	nrnodes=2*max(nsample,ntest)+1,
     &	mimp=imp*(mdim-1)+1,
     &	nimp=impn*(nsample-1)+1,
     &	nprot0=max(nprot,1),
     &	nscale0=max(nscale,1),
     &	ntest0=max(ntest,1),
     &	nprdata=min(nprox,1),
     &	near=nprdata*(nsample-1)+1,
     &	nearts=nprdata*(ntest0-1)+1)
c
c	-------------------------------------------------------
c	DIMENSIONING OF ARRAYS
c
	real x(mdim,nsample),xts(mdim,ntest0),fill(mdim,nclass),
     &	tgini(mdim),tsgini(mdim),zt(mdim),avgini(mdim),rfill(mdim),
     &  xxts(mdim,nclass*ntest0),yg(mdim,nprot0,nclass),
     &	votes(mdim,jbt),effect(mdim,mdim),teffect(mdim,mdim),
     &	hist(0:mdim,mdim),g(mdim),
     &	dgini(nrnodes),xbestsplit(nrnodes),tnodewt(nrnodes),
     &	tw(nrnodes),tn(nrnodes),
     &	v(nsample),win(nsample),
     &	q(nclass,nsample),tk(0:nprot0,nsample),devout(nclass),
     &	classwt(nclass),wr(nclass),tmissts(nclass), ff(nclass),
     &	tmiss(nclass),tclasspop(nclass),wl(nclass),rmedout(nclass),
     &	tclasscat(nclass,maxcat),qts(nclass,ntest0),
     &	qxxts(nclass,nclass*ntest0),classpop(nclass,nrnodes),
     &	signif(mimp),zscore(mimp),sqsd(mimp),avimp(mimp),
     &	qimp(nclass,nimp,mimp),
     &	tout(near),outtr(near),
     &	xc(maxcat),dn(maxcat),cp(maxcat),cm(maxcat),
     &	votecat(maxcat),fp(nprot0,nclass),
     &	outts(nearts),touts(nearts)
c
    	integer cat(mdim),mbax(mdim),iv(mdim),msm(mdim),muse(mdim),
     &	kf(mdim,nclass),irnk(mdim,jbt),missing(mdim,near),
     &	a(mdim,nsample),at(mdim,nsample),b(mdim,nsample),
     &	cl(nsample),out(nsample),nodextr(nsample),nodexvr(nsample),
     &	jin(nsample),joob(nsample),pjoob(nsample),
     &	jvr(nsample),jtr(nsample),jest(nsample),
     &	isort(nsample),ncass(nsample0),
     &	ta(nsample),ncase(nsample),idmove(nsample),
     &	jests(ntest0),jints(ntest0),jts(ntest0),
     &	nodexts(ntest0),clts(ntest0),imax(ntest0),
     &	bestsplitnext(nrnodes),bestvar(nrnodes),bestsplit(nrnodes),
     &	nodestatus(nrnodes),nodepop(nrnodes),nodestart(nrnodes),
     &	nodeclass(nrnodes),
     &	parent(nrnodes),treemap(2,nrnodes),
     &	ncts(nclass),nc(nclass),mtab(nclass,nclass),jl(nclass),
     &	jxxts(nclass*ntest0),nodexxts(nclass*ntest0),
     &	ndpointxxs(nclass*ntest0,jbt),
     &	nrcat(maxcat),krcat(maxcat),kcat(maxcat),ncatsplit(maxcat),
     &	nbestcat(maxcat,nrnodes),
     &  kpx(near),kpr(near),jpr(near),inuse(near),key(near),
     &	nodex(near),kend(near),ncp(near),
     &	ndpoint(near,jbt),jprox(near,jbt),
     &	npcase(near,jbt),kprox(near,jbt),loz(near,nrnn),
     &	iout(nearts),kendts(nearts),ncpts(nearts),
     &	ndpoints(nearts,jbt),lozts(nearts,nrnn),
     &	lcase(0:jbt),kprot(nprot0,nclass)
c
	character*500 text
c
c	-------------------------------------------------------
c    	SCALAR DECLARATIONS
c
	real right, rightimp, errtr, errts, tavg, er, zz, rrand 
c
	integer mtry, last, nr, n, 
     &	m, mdimt, k, j, nrep, i, m1, kv, jp, ncat, nn,
     &	mt, jb, nuse, ndbigtree, nout, jj, nts, mr,
     &	n0, n1, n2, n3, n4, n5, n6, n7, n8, n9
c
c	-------------------------------------------------------
c    	USED IN PROXIMITY AND SCALING
c
	double precision prox(near,nrnn),y(near),u(near),dl(nscale0),
     &	xsc(near,nscale0),red(near),ee(near),ev(near,nscale0),
     &	proxts(nearts,nrnn)
c
c	-------------------------------------------------------
c	READ OLD TREE STRUCTURE AND/OR PARAMETERS
c
 	if (irunrf.eq.1)
     &		open(1,file='savedforest',status='old')
	if (ireadpar.eq.1)
     &		open(2,file='savedparams',status='old')
	if (ireadfill.eq.1)
     &		open(3, file='savedmissfill',status='old')
 	if (ireadprox.eq.1)
     &		open(4, file='savedprox', status='old')
c
c	-------------------------------------------------------
c	NAME OUTPUT FILES FOR SAVING THE FOREST STRUCTURE
c
 	if (isaverf.eq.1)
     &		open(1, file='savedforest',status='new')
	if (isavepar.eq.1)
     &		open(2, file='savedparams',status='new')
	if (isavefill.eq.1)
     &		open(3, file='savedmissfill',status='new')
 	if (isaveprox.eq.1)
     &		open(4, file='savedprox', status='new')
c
c	-------------------------------------------------------
c	NAME OUTPUT FILES TO SAVE DATA FROM CURRENT RUN
c
 	if (idataout.ge.1)
     &		open(7, file='save-data-from-run',status='new')
	if (impfastout.eq.1)
     &		open(8,file='save-impfast',status='new')
 	if (impout.eq.1)
     &		open(9,file='save-importance-data',status='new')
 	if (impnout.eq.1)
     &		open(10,file='save-caseimp-data',status='new')
	if (interout.eq.1)
     &		open(11,file='save-pairwise-effects',status='new')
	if (iprotout.eq.1)
     &		open(12,file='save-protos',status='new')
 	if (iproxout.ge.1)
     &		open(13, file='save-run-proximities',status='new')
	if (iscaleout.eq.1)
     &		open(14, file='save-scale',status='new')
	if (ioutlierout.ge.1)
     &		open(15, file='save-outliers',status='new')
c
c	-------------------------------------------------------
c	READ IN DATA--SEE MANUAL FOR FORMAT
c
	open(16, file='satimage.tra', status='old')
	do n=1,nsample0
		read(16,*) (x(m,n),m=1,36), cl(n)
	end do
	close(16)
	if(ntest.gt.1) then
		open(17, file='satimage.tes', status='old')
		do n=1,ntest0
			read(17,*) (xts(m,n),m=1,36),clts(n)
		end do
		close(17)
	end if
c
c	-------------------------------------------------------
c	SELECT SUBSET OF VARIABLES TO USE
c
	if(mselect.eq.0) then
		mdimt=mdim
		do k=1,mdim
			msm(k)=k
		enddo
	endif
c	if (mselect.eq.1) then
c		fill in which variables 
c		mdimt=
c		msm(1)=
c		---
c		msm(mdimt)=
c	endif
c
c	-------------------------------------------------------
c	SET CATEGORICAL VALUES
c
	do m=1,mdim
		cat(m)=1
	enddo
c	fill in cat(m) for all variables m for which cat(m)>1
c	do m=1,mdim
c		cat(m)=4
c	enddo
c
c	-------------------------------------------------------
c	SET CLASS WEIGHTS
c
	do j=1,nclass
		classwt(j)=1
	enddo
c	if(jclasswt.eq.1) then
c		fill in classwt(j) for each j:
c		classwt(1)=1.
c		classwt(2)=10.
c	endif
c
c	=======================================================
c	**************  END OF USER INPUT  ********************
c	=======================================================
c
c	-------------------------------------------------------
c	SPECIAL CASES - QUERY PARAMETERS OR USE SAVED TREE
c		
c	query parameters
	if(ireadpar.eq.1) goto 888
c		
c	use saved tree
	if(irunrf.eq.1) goto 999
c
c	-------------------------------------------------------
c	ERROR CHECKING
c
	if(labelts.ne.0.and.labelts.ne.1) then
		write(*,*) "error labelts",labelts
		stop
	endif
	if(labeltr.ne.0.and.labeltr.ne.1) then
		write(*,*) "error labeltr",labeltr
		stop
	endif
	if(labeltr.eq.0.and.nclass.ne.2) then
		write(*,*) "error, nclass should be 2 if labeltr=0"
		stop
	endif
	if(lookcls.ne.0.and.lookcls.ne.1) then
		write(*,*) "error lookcls",lookcls
		stop
	endif
	if(jclasswt.ne.0.and.jclasswt.ne.1) then
		write(*,*) "error jclasswt",jclasswt
		stop
	endif
	if(mselect.ne.0.and.mselect.ne.1) then
		write(*,*) "error mselect",mselect
		stop
	endif
	if(mdim2nd.lt.0.or.mdim2nd.gt.mdim) then
		write(*,*) "error mdim2nd",mdim2nd
		stop
	endif
	if(imp.ne.0.and.imp.ne.1) then
		write(*,*) "error imp",imp
		stop
	endif
	if(impn.ne.0.and.impn.ne.1) then
		write(*,*) "error impn",impn
		stop
	endif
	if(interact.ne.0.and.interact.ne.1) then
		write(*,*) "error interact",interact
		stop
	endif
	if(impfast.ne.0.and.impfast.ne.1) then
		write(*,*) "error impfast",impfast
		stop
	endif
	if(nprox.lt.0.or.nprox.gt.2) then
		write(*,*) "error nprox",nprox
		stop
	endif
	if(nrnn.lt.0.or.nrnn.gt.nsample) then
		write(*,*) "error - nrnn",nrnn
		stop
	endif
	if(noutlier.lt.0.or.noutlier.gt.2) then
		write(*,*) "error - noutlier",noutlier
		stop
	endif
	if(nscale.lt.0.or.nscale.gt.mdim) then
		write(*,*) "error - nscale",nscale
		stop
	endif
	if(nprot.lt.0.or.nprot.gt.nsample) then
		write(*,*) "error - nprot",nprot
		stop
	endif
	if(mfixall.lt.0.or.mfixall.gt.2) then
		write(*,*) "error - mfixall",mfixall
		stop
	endif

	if(isaverf.ne.0.and.isaverf.ne.1) then
		write(*,*) "error isaverf",isaverf
		stop
	endif
	if(irunrf.ne.0.and.irunrf.ne.1) then
		write(*,*) "error irunrf",irunrf
		stop
	endif
	if(isavepar.ne.0.and.isavepar.ne.1) then
		write(*,*) "error isavepar",isavepar
		stop
	endif
	if(ireadpar.ne.0.and.ireadpar.ne.1) then
		write(*,*) "error ireadpar",ireadpar
		stop
	endif
	if(isavefill.ne.0.and.isavefill.ne.1) then
		write(*,*) "error isavefill",isavefill
		stop
	endif
	if(ireadfill.ne.0.and.ireadfill.ne.1) then
		write(*,*) "error ireadfill",ireadfill
		stop
	endif
	if(isaveprox.ne.0.and.isaveprox.ne.1) then
		write(*,*) "error isaveprox",isaveprox
		stop
	endif
	if(ireadprox.ne.0.and.ireadprox.ne.1) then
		write(*,*) "error ireadprox",ireadprox
		stop
	endif
	if(isumout.lt.0.or.isumout.gt.1) then
		write(*,*) "error isumout",isumout
		stop
	endif
	if(idataout.lt.0.or.idataout.gt.2) then
		write(*,*) "error idataout",idataout
		stop
	endif
	if(impfastout.lt.0.or.impfastout.gt.1) then
		write(*,*) "error impfastout",impfastout
		stop
	endif
	if(impout.lt.0.or.impout.gt.2) then
		write(*,*) "error impout",impout
		stop
	endif
	if(interout.lt.0.or.interout.gt.2) then
		write(*,*) "error interout",interout
		stop
	endif
	if(iprotout.lt.0.or.iprotout.gt.2) then
		write(*,*) "error iprotout",iprotout
		stop
	endif
	if(iproxout.lt.0.or.iproxout.gt.2) then
		write(*,*) "error iproxout",iproxout
		stop
	endif
	if(iscaleout.lt.0.or.iscaleout.gt.1) then
		write(*,*) "error iscaleout",iscaleout
		stop
	endif
	if(ioutlierout.lt.0.or.ioutlierout.gt.2) then
		write(*,*) "error ioutlierout",ioutlierout
		stop
	endif
	if(mfixtr.gt.0.and.mfixall.gt.0) then
		write(*,*) "error - can't do both careful and rough fix"
		write(*,*) "mfixtr=",mfixtr
		write(*,*) "mfixall=",mfixall
		stop
	endif
	if(mfixtr.gt.0.and.nprox.eq.0) then
		write(*,*) "error - mfixtr>0 and nprox=0"
		stop
	endif
	if(noutlier.gt.0.and.nprox.eq.0) then
		write(*,*) "error - noutlier>0 and nprox=0"
		stop
	endif
	if(nscale.gt.0.and.nprox.eq.0) then
		write(*,*) "error - nscale>0 and nprox=0"
		stop
	endif
	if(nprot.gt.0.and.nprox.eq.0) then
		write(*,*) "error - nprot>0 and nprox=0"
		stop
	endif
	if(mdim2nd.gt.0.and.imp.eq.0) then
		write(*,*) "error - mdim2nd>0 and imp=0"
		stop
	endif
	if(impn.gt.0.and.imp.eq.0) then
		write(*,*) "error - impn>0 and imp=0"
		stop
	endif
	do m=1,mdim
		if(cat(m).gt.maxcat) then
			write(*,*) "error in cat",m,cat(m)
			stop
		endif
	enddo
	if(labeltr.eq.1) then
		do n=1,nsample
			if(cl(n).lt.1.or.cl(n).gt.nclass) then 
				write(*,*) 'error in class label', n,cl(n)
				stop
			endif
		enddo
	endif
c
c	-------------------------------------------------------
c	INITIALIZATION
c
c	mtry will change if mdim2nd>0
	mtry=mtry0
c
c	last = 1 if this is the final number of variables
	if(mdim2nd.gt.0) then
		last=0
	else
		last=1
	endif
c
c	assign node weights for equal weighting case (use getweights 
c	for unequal weighting case)
	if(jclasswt.eq.0) then
		do k=1,nrnodes
			tnodewt(k)=1
		enddo
	endif
c
c	set seed
	do nr=1,111
		zz=rrand(1)
	enddo
c
c	-------------------------------------------------------
c	SET UP DATA TO ADD A CLASS
c
	if(labeltr.eq.0) then
		call createclass(x,cl,nsample0,nsample,mdim)
	endif
c
c	-------------------------------------------------------
c	COUNT CLASS POPULATIONS
c
	call zerv(nc,nclass)
	do n=1,nsample
		nc(cl(n))=nc(cl(n))+1
	enddo
	if(ntest.gt.0.and.labelts.eq.1) then
		call zerv(ncts,nclass)
		do n=1,ntest0
			ncts(clts(n))=ncts(clts(n))+1
		enddo
	endif
c
c	-------------------------------------------------------
c	DO PRELIMINARY MISSING DATA
c
	if(mfixtr.gt.0) then
		nrep=1
		call roughfix(x,v,mdim,nsample,near,
     &		missing,cat,code,nrcat,maxcat,rfill)
	else
		nrep=0
		if(mfixall.ge.1) then
c			fast fix on the training data
			if(ireadfill.eq.1) then
				read(3,*)((fill(m,j),m=1,mdim),j=1,nclass)
			else
				call fixfast(x,v,ncase,mdim,nsample,cat,code,
     &				nrcat,maxcat,fill,cl,nclass,kf,ff,jl,krcat)
			endif
			call xfill(x,nsample,mdim,fill,cl,nclass,code)
		endif
		if(mfixall.eq.2.and.ntest.gt.0) then
c			fast fix on the test data - different method if labeled
			if(labelts.eq.0) then
				call fixtestunlab(xxts,xts,mdim,nclass,ntest,fill,code)
			else
				call fixtestlab(xts,mdim,nclass,ntest,fill,code,clts)
			endif
		endif
		if(isavefill.eq.1) then
			write(3,*)((fill(m,j),m=1,mdim),j=1,nclass)
		endif
	endif
c
	if(ireadprox.eq.1) then
c		read in saved proximities
		do jb=1,jbt
			do n=1,near
				read(4,*) jprox(n,jb),kprox(n,jb)
			enddo
		enddo
	endif
c
c	-------------------------------------------------------
c	TOP OF LOOP FOR ITERATIVE MISSING VALUE FILL OR SUBSET SELECTION
c	i.e. mfixtr>0 or mdim2nd>0
c
13	continue 
c
c	-------------------------------------------------------
c	INITIALIZE FOR RUN
c
	call zermr(q,nclass,nsample)
	call zerv(out,nsample)
	call zervr(avgini,mdim)
	call zervr(tsgini,mdim)
	call zermd(prox,near,nrnn)
	if(imp.gt.0) then
		call zervr(sqsd,mimp)
		call zervr(avimp,mimp)
	endif
	call zermr(votes,mdim,jbt)
	if(ntest.gt.0) then
 		call zermr(qts,nclass,ntest)
		call zerv(jints,ntest)
		call zermr(qxxts,nclass,nclass*ntest)
	endif
	if(impn.eq.1) then
		do mt=1,mdimt
			m=msm(mt)
			do n=1,nsample
				do k=1,nclass
					qimp(k,n,m)=0
				enddo
			enddo
		enddo
	endif
	call makea(x,mdim,nsample,cat,isort,v,at,b,mbax,mdimt,msm)
c
c	=======================================================
c	**************  BEGIN MAIN LOOP  **********************
c	=======================================================
c
	do jb=1,jbt
c
c	-------------------------------------------------------
c	INITIALIZE
c
	call zerv(nodestatus,nrnodes)
	call zerm(treemap,2,nrnodes)
	call zervr(xbestsplit,nrnodes)
	call zerv(nodeclass,nrnodes)
	call zerv(bestvar,nrnodes)
	call zerv(jin,nsample)
	call zervr(tclasspop,nclass)
	call zervr(win,nsample)
c
c	-------------------------------------------------------
c	TAKE A BOOTSTRAP SAMPLE 
c
	do n=1,nsample
		k=int(rrand(1)*nsample)+1
		win(k)=win(k)+classwt(cl(k))
		jin(k)=jin(k)+1
		tclasspop(cl(k))=tclasspop(cl(k))+classwt(cl(k))
	enddo
c
c	-------------------------------------------------------
c	PREPARE TO BUILD TREE
c
	call eqm(a,at,mdim,nsample,mdimt,msm)
	call moda(a,nuse,nsample,mdim,cat,maxcat,ncase,jin,mdimt,msm)
c
c	-------------------------------------------------------
c	MAIN TREE-BUILDING 
c
	call buildtree(a,b,cl,cat,mdim,nsample,nclass,treemap,
     &  bestvar,bestsplit,bestsplitnext,dgini,nodestatus,nodepop,
     &	nodestart,classpop,tclasspop,tclasscat,ta,nrnodes,
     &	idmove,ndsize,ncase,parent,mtry,nodeclass,ndbigtree,
     &	win,wr,wl,nuse,kcat,ncatsplit,xc,dn,cp,cm,maxcat,
     &	nbestcat,msm,mdimt)
c
c	-------------------------------------------------------
c	SPLIT X
c
	call xtranslate(x,mdim,nrnodes,nsample,bestvar,bestsplit,
     &  bestsplitnext,xbestsplit,nodestatus,cat,ndbigtree)
c
c	-------------------------------------------------------
c	ASSIGN CLASSWEIGHTS TO NODES
c
	if(jclasswt.eq.1) then
		call getweights(x,nsample,mdim,treemap,nodestatus,
     &		xbestsplit,bestvar,nrnodes,ndbigtree,
     &		cat,maxcat,nbestcat,jin,win,tw,tn,tnodewt)
	endif
c
c	-------------------------------------------------------
c	GET OUT-OF-BAG ESTIMATES
c
	call testreebag(x,nsample,mdim,treemap,nodestatus,
     &  xbestsplit,bestvar,nodeclass,nrnodes,ndbigtree,
     &	cat,jtr,nodextr,maxcat,nbestcat,dgini,tgini)
c
	do j=1,mdimt
		m = msm(j)
		votes(m,jb)=tgini(m)
		avgini(m)=avgini(m)+tgini(m)
	enddo
c
	nout=0
	right=0
	do n=1,nsample
		if(jin(n).eq.0) then
c			this case is out-of-bag
			q(jtr(n),n)=q(jtr(n),n)+tnodewt(nodextr(n))
c			update count of correct oob classifications
			if(jtr(n).eq.cl(n)) right=right+tnodewt(nodextr(n))
c			out(n) = number of TREES for which the nth obs has been out-of-bag
			out(n) = out(n)+1
c			nout = total number of obs that are out-of-bag for THIS tree
			nout=nout+1
			joob(nout) = n
		endif
	enddo
c
c	-------------------------------------------------------
c	GET TEST SET ERROR ESTIMATES
c
	if(ntest.gt.0) then
		if(labelts.eq.1.or.mfixall.lt.2) then
c			labelled test set or no missing value fills - use xts
			call testreebag(xts,ntest,mdim,treemap,nodestatus,
     & 			xbestsplit,bestvar,nodeclass,nrnodes,ndbigtree,
     &			cat,jts,nodexts,maxcat,nbestcat,dgini,tsgini)
c
			do n=1,ntest0
				qts(jts(n),n)=qts(jts(n),n)+tnodewt(nodexts(n))
			enddo
		endif
		if(labelts.eq.0.and.mfixall.eq.2) then
c			unlabelled test set with missing values - use xxts
			nts=nclass*ntest0
			call testreebag(xxts,nts,mdim,treemap,nodestatus,
     & 			xbestsplit,bestvar,nodeclass,nrnodes,ndbigtree,
     &			cat,jxxts,nodexxts,maxcat,nbestcat,dgini,tsgini)
c
			do n=1,nclass*ntest0
				k=(n/nclass)+1
				qxxts(jxxts(n),n)=qxxts(jxxts(n),n)+tnodewt(nodexxts(k))
			enddo
		endif
	endif
c
c	-------------------------------------------------------
c	PRE-PROX 
c
	if(nprox.ge.1) then
		do n=1,near
			ndpoint(n,jb)=nodextr(n)
			nodex(n)=nodextr(n)
		enddo
		call preprox(nodex,near,ncass,jpr,kpx)
		do n=1,near
			if(ireadprox.eq.0) then
				jprox(n,jb)=jpr(n)
				kprox(n,jb)=kpx(n)
			endif
			npcase(n,jb)=ncass(n)
		enddo
		do n=1,nearts
			ndpoints(n,jb)=nodexts(n)
		enddo
	endif
c
c	-------------------------------------------------------
c	VARIABLE IMPORTANCE
c	(ONLY IF WE ARE DONE WITH IMPUTATION)
c
	if(imp.eq.1.and.nrep.eq.mfixtr) then
		call zerv(iv,mdim)
		do jj=1,ndbigtree
c			iv(j)=1 if variable j was used to split on
			if(nodestatus(jj).ne.-1) iv(bestvar(jj))=1
		enddo
		do k=1,mdimt
			mr=msm(k)
c			choose only those that used a split on variable mr 
c
			if(iv(mr).eq.1) then
				call testreeimp(x,nsample,mdim,joob,pjoob,nout,mr,
     &				treemap,nodestatus,xbestsplit,bestvar,nodeclass,nrnodes,
     &				ndbigtree,cat,jvr,nodexvr,maxcat,nbestcat)
c
				rightimp=0
				do n=1,nout
c					the nth out-of-bag case is the nnth original case
					nn = joob(n)
					if(impn.eq.1) then
						qimp(jvr(n),nn,mr)=qimp(jvr(n),nn,mr)+
     &						tnodewt(nodexvr(n))
					endif
					if(jvr(n).eq.cl(nn)) rightimp=rightimp+tnodewt(nodexvr(n))
				enddo
				avimp(mr)=avimp(mr)+(right-rightimp)/nout
				sqsd(mr)=sqsd(mr)+((right-rightimp)**2)/nout
			else
				if(impn.eq.1) then
					do n=1,nsample
						if(jin(n).eq.0) then
							qimp(jtr(n),n,mr)=qimp(jtr(n),n,mr)+
     &							tnodewt(nodextr(n))
						endif
					enddo
				endif
			endif
		enddo !k
	endif
c
c	-------------------------------------------------------
c	GIVE RUNNING OUTPUT
c
	if(lookcls.eq.1.and.jb.eq.1) then
		write(*,*) 'class counts-training data'
		write(*,*) (nc(j), j=1,nclass)
		if(ntest.gt.0.and.labelts.eq.1) then
			print*
			write(*,*) 'class counts-test data'
			write(*,*) (ncts(j), j=1,nclass)
		endif
	endif
	if(mod(jb,look).eq.0.or.jb.eq.jbt) then
		call comperrtr(q,cl,nsample,nclass,errtr,
     &  	tmiss,nc,jest,out)
		if(lookcls.eq.1) then
			write(*,'(i8,100f10.2)') 
     &			jb,100*errtr,(100*tmiss(j), j=1,nclass)
		else
			write(*,'(i8,2f10.2)') jb,100*errtr
		endif
		if(ntest.gt.0.and.labelts.eq.1) then
			call comperrts(qts,clts,ntest,nclass,errts,
     & 			tmissts,ncts,jests)
			if(lookcls.eq.1) then
				write(*,'(i8,20f10.2)') jb,100*errts,
     &				(100*tmissts(j),j=1,nclass)
			else
				write(*,'(i8,2f10.2)') jb,100*errts
			endif
			print *
		endif
	endif
c
c	-------------------------------------------------------
c	SEND SAVETREE DATA TO FILE (IF THIS IS THE FINAL FOREST)
c
	if(isaverf.eq.1.and.last.eq.1) then
		if(jb.eq.1) write(1,*) (cat(m),m=1,mdim)
		write(1,*) ndbigtree
   		do n=1,ndbigtree
			write(1,*) n,nodestatus(n),bestvar(n),treemap(1,n),
     &			treemap(2,n),nodeclass(n),xbestsplit(n),tnodewt(n),
     &			(nbestcat(k,n),k=1,maxcat)
		enddo
	endif
c
c	=======================================================
c	**************  END MAIN LOOP  ************************
c	=======================================================
	enddo !jb
c
c	-------------------------------------------------------
c	FIND PROXIMITIES, FILL IN MISSING VALUES AND ITERATE
c
	if(nrep.lt.mfixtr) then
c		not done (note: if mfixtr=0, nrep=0)
		nrep=nrep+1
		write(*,*) 'nrep', nrep
c		compute proximities between each obs in the training 
c		set and its nrnn closest neighbors
		call comprox(kprox,ndpoint,jprox,npcase,kpr,
     &		lcase,loz,prox,near,jbt,nrnn,kend,key,inuse)
c		use proximities to impute missing values
		call impute(x,prox,near,mdim,nclass,
     &		maxcat,votecat,cat,nrnn,loz,cl,kend,
     &		missing)
c		go to the top of the main loop
		goto 13
	endif
c
c	-------------------------------------------------------
c	COMPUTE IMPORTANCES, SELECT SUBSET AND LOOP BACK
c
	if(imp.eq.1) then
		call finishimp(mdim,sqsd,avimp,signif,
     &		zscore,jbt,mdimt,msm)
		do m=1,mdimt
c			zt(m) is the negative z-score for variable msm(m)
			zt(m)=-zscore(msm(m))
			muse(m)=msm(m)
		enddo
c		sort zt from smallest to largest, and sort muse accordingly
		call quicksort(zt,muse,1,mdimt,mdim)
c		muse(m) refers to the variable that has the mth-smallest zt
c		select the mdim2nd most important variables and iterate
		if(last.eq.0.and.mdim2nd.ge.1) then
			mdimt=mdim2nd
			do m=1,mdimt
				msm(m)=muse(m)
			enddo
			last=1 !done variable selection
			mtry=nint(sqrt(real(mdimt)))
			goto 13
	   	endif
	endif
c
c	-------------------------------------------------------
c	END OF ITERATIONS - NOW ENDGAME
c	-------------------------------------------------------
c
c	-------------------------------------------------------
c	NORMALIZE VOTES
c
	do j=1,nclass
		do n=1,nsample
			if(q(j,n).gt.0.0.and.out(n).gt.0) q(j,n)=q(j,n)/real(out(n))
		enddo
		if(ntest.gt.0) then
			do n=1,ntest0
				qts(j,n)=qts(j,n)/jbt
			enddo
		endif
	enddo
c
c	-------------------------------------------------------
c	COMPUTE PROXIMITIES AND SEND TO FILE
c
	if(nprox.ge.1) then
c		compute proximities between each obs in the training 
c		set and its nrnn closest neighbors
		call comprox(kprox,ndpoint,jprox,npcase,kpr,
     &		lcase,loz,prox,near,jbt,nrnn,kend,key,inuse)
		if(iproxout.ge.1) then 
			do n=1,near
				write(13,'(i5,500(i5,f10.3))') n,(loz(n,k),
     &				prox(n,k), k=1,nrnn)
			enddo
		endif
c		save proximities for use in future tree
		if(isaveprox.eq.1) then
			do jb=1,jbt
				do n=1,near
					write(4,*) jprox(n,jb),kprox(n,jb)
				enddo
			enddo
		endif
    	endif
	if(nprox.eq.2.and.ntest.gt.0) then
c		compute proximities between each obs in the test 
c		set and its nrnn closest neighbors in the training set
		call comproxts(kprox,jprox,npcase,kpr,lcase,near,
     &		jbt,nrnn,nearts,proxts,lozts,ndpoints,kendts,key,inuse)
		if(iproxout.eq.2) then 
			do n=1,nearts
				write(13,'(i5,500(i5,f10.3))') n,(lozts(n,k),
     &				proxts(n,k), k=1,nrnn)
			enddo
		endif
	endif
	close(13)
c
c	-------------------------------------------------------
c	COMPUTE SCALING COORDINATES AND SEND TO FILE
c
	if(nprox.eq.1.and.nscale.gt.0) then
		call myscale(loz,prox,xsc,y,u,near,nscale,red,nrnn,
     &		ee,ev,dl,kend)
		if(iscaleout.eq.1)then
			do n=1,nscale0
				write(14,'(i5,f10.3)') n,dl(n)
			enddo
			do n=1,near
				write(14,'(3i5,15f10.3)') n,cl(n),jest(n),
     &				(xsc(n,k),k=1,nscale0)
			enddo
			close(14)
		endif
	endif
c
c	-------------------------------------------------------
c	COMPUTE CASEWISE VARIABLE IMPORTANCE (FOR GRAPHICS)
c	AND SEND TO FILE
c
	if(impn.eq.1.and.impnout.eq.1) then
		do m1=1,mdimt
			mr=msm(m1)
			do n=1,nsample
				do j=1,nclass
					qimp(j,n,mr)=q(j,n)-(qimp(j,n,mr)/out(n))
				enddo
			enddo
		enddo
		do j=1,nclass
			do n=1,nsample
				write(10,*)(qimp(j,n,msm(m1)),m1=1,mdimt)
			enddo
		enddo
	endif
c
c	-------------------------------------------------------
c	SEND IMPORTANCES TO FILE OR SCREEN
c
	if(imp.eq.1) then
		do k=1,min0(mdimt,25)
			m=muse(k)
			if(impout.eq.1) write(9,'(i5,10f10.3)')m,100*avimp(m), 
     &			zscore(m),signif(m)
			if(impout.eq.2) write(*,'(i5,10f10.3)')m,100*avimp(m), 
     &			zscore(m),signif(m)
		enddo
		close(9)
	endif
c
c	-------------------------------------------------------
c	COMPUTE INTERACTIONS AND SEND TO FILE OR SCREEN
c
	if(interact.eq.1) then
		call compinteract(votes,effect,msm,mdim,mdimt,
     &		jbt,g,iv,irnk,hist,teffect)
		if(interout.eq.1) then
			write(11,*)'CODE'
			do i=1,mdimt
				write(11,*) i,msm(i)
			enddo
			print*
			do i=1,mdimt
				m=msm(i)
				effect(m,m)=0
				write(11,'(40i5)') i,
     &				(nint(effect(m,msm(j))), j=1,mdimt)
			enddo
			close(11)
		endif
		if(interout.eq.2) then
			write(*,*)'CODE'
			do i=1,mdimt
				write(*,*) i,msm(i)
			enddo
			print*
			do i=1,mdimt
				m=msm(i)
				effect(m,m)=0
				write(*,'(20i5)') i,
     &				(nint(effect(m,msm(j))), j=1,mdimt)
			enddo
		endif
	endif
c
c	-------------------------------------------------------
c	COMPUTE FASTIMP AND SEND TO FILE
c
	if(impfast.eq.1.and.impfastout.eq.1) then
		write(*,*) "here"
		tavg=0
		do k=1,mdimt
			m=msm(k)
			tavg=tavg+avgini(m)
		enddo
		tavg=tavg/mdimt
		do k=1,mdimt
			m=msm(k)
			write(8,*) m,avgini(m)/tavg
		enddo
		close(8)
	endif
c
c	-------------------------------------------------------
c	FINISH MFIXALL=2 AND LABELTS=0
c
	if(labelts.eq.0.and.mfixall.eq.2.and.ntest.gt.0) then
		call reducemult(wr,nclass,ntest,qxxts,jests,qts,
     &		xxts,xts,mdim,imax)
	endif
c
c	-------------------------------------------------------
c	COMPUTE PROTOTYPES AND SEND TO FILE
c
	if(nprox.ge.1.and.nprot.gt.0) then
		call compproto(prox,jest,fp,loz,nsample,nprot,tk,nrnn,
     &		nclass,kprot,kend)
c
		do kv=1,mdimt
			m=msm(kv)
			if(cat(m).eq.1) then
				do jp=1,nclass
					do k=1,nprot0
						yg(m,k,jp)=real(b(m,kprot(k,jp)))/mbax(m)
					enddo
				enddo
			else
				ncat=cat(m)
				do jp=1,nclass
					do k=1,nprot0
						yg(m,k,jp)=(x(m,kprot(k,jp))-1)/(ncat-1)
					enddo
				enddo
			endif
		enddo
		if(iprotout.eq.1)write(12,'(i5,50f7.2)') 0,
     &		((fp(k,jp), k=1,nprot0),jp=1,nclass)
		if(iprotout.eq.2)write(*,'(i5,50f7.2)') 0,
     &		((fp(k,jp), k=1,nprot0),jp=1,nclass)
		do kv=1,mdimt
			if(iprotout.eq.1) write(12,'(i5,50f7.2)')kv,
     &			((yg(msm(kv),k,jp), k=1,nprot0),jp=1,nclass)
			if(iprotout.eq.2) write(*,'(i5,50f7.2)')kv,
     &			((yg(msm(kv),k,jp), k=1,nprot0),jp=1,nclass)
    		enddo
		close (12)
	endif
c
c	-------------------------------------------------------
c	COMPUTE OUTLIER MEASURE AND SEND TO FILE
c
	if(noutlier.ge.1) then
		call locateout(cl,near,nclass,
     &		labeltr,outtr,tout,isort,devout,rmedout,ncp,
     &		ndpoint,jprox,npcase,kpr,kprox,jbt)
		do n=1,near
			if(ioutlierout.ge.1) write(15,'(2i5,f10.3)') n,cl(n),
     &			amax1(outtr(n),0.0)
		enddo
	endif
	if(noutlier.eq.2) then
		call locateouts(clts,nearts,labelts,
     &		outts,jests,ncpts,ndpoints,jprox,npcase,kpr,kprox,
     &		near,jbt,touts,iout)
		do n=1,nearts
			if(ioutlierout.eq.2) write(15,'(3i5,f10.3)') n,jests(n),
     &			amax1(outts(n),0.0)
		enddo
		close(15)
	endif
c
c	-------------------------------------------------------
c	SUMMARY OUTPUT
c
	if (isumout.eq.1) then
		if(ntest.eq.0.or.labelts.eq.0)then
			write(*,*) 'final error rate %    ',100*errtr
		else
			write(*,*) 'final error rate %    ',100*errtr
			write(*,*) 'final error test %    ',100*errts
		endif
		print *
		call zerm(mtab,nclass,nclass)
		do n=1,nsample
			if(jest(n).gt.0) mtab(cl(n),jest(n))=mtab(cl(n),jest(n))+1
		enddo
		write(*,*) '          true class '
		print *
		write(*,'(20i6)')  (i,i=1,nclass)
		print *
		do j=1,nclass
			write(*,'(20i6)')  j,(mtab(i,j), i=1,nclass)
		enddo
		print *
	endif
c
c	-------------------------------------------------------
c	SEND INFO ON TRAINING AND/OR TEST SET DATA TO FILE 
c
	if(idataout.ge.1) then
		do n=1,nsample
			write(7,'(3i5,1000f10.3)') n,cl(n),jest(n),
     & 			(q(j,n), j=1,nclass),(x(m,n),m=1,mdim)
		enddo
	endif
	if(idataout.eq.2.and.ntest.gt.0) then
		if(labelts.eq.1) then
			do n=1,ntest0
				write(7,'(3i5,1000f10.3)') n,clts(n),jests(n),
     & 				(qts(j,n), j=1,nclass),(xts(m,n),m=1,mdim)
			enddo
		else
			do n=1,ntest0
				write(7,'(3i5,1000f10.3)') n,jests(n),
     & 				(qts(j,n), j=1,nclass),(xts(m,n),m=1,mdim)
			enddo
		endif
	endif
	close(7)
c
c	-------------------------------------------------------
c	SEND RUN PARAMETERS TO FILE 
c
	if(isavepar.eq.1) then
		write(2,*) nsample0,mdim,maxcat,nclass,jbt,
     &		jclasswt,mfixall,mfixtr,code,nrnodes,
     &		100*errtr
c
c		type in comments up to 500 characters long
c		between the ' ' in the line below.
c
		write(2,*) 'this is a test run to verify that my 
     &		descriptive output works.'
		close (2)
	endif
c
c	-------------------------------------------------------
c	END OF USUAL RUN
c
	stop
c
c	-------------------------------------------------------
c	READ RUN PARAMETERS AND PRINT TO SCREEN
c
888	if(ireadpar.eq.1) then
		read(2,*) n0,n1,n2,n3,n4,n5,n6,n7,n8,n9,er
		write(*,*) 'parameters'
		write(*,*) 'nsample=' ,n0
		write(*,*) 'mdim =  ' ,n1
		write(*,*) 'maxcat = ',n2
		write(*,*) 'nclass = ',n3
		write(*,*) 'jbt =   ' ,n4
		write(*,*) 'jclasswt=',n5
		write(*,*) 'mfixall=' ,n6
		write(*,*) 'mfixtr='  ,n7
		write(*,*) 'code='    ,n8
		write(*,*) 'nrnodes=' ,n9
		print *
		write(*,*) 'out-of-bag error =', er,'%'
		print *
		read(2,'(500a)') text
		write(*,*) text
c
		stop
c
	endif
c
c	-------------------------------------------------------
c	RERUN OLD RANDOM FOREST
c
999	if(irunrf.eq.1.and.ntest.gt.0) then
		call runforest(mdim,ntest,nclass,maxcat,nrnodes,
     &		labelts,jbt,clts,xts,xbestsplit,qts,treemap,nbestcat,
     &		nodestatus,cat,nodeclass,jts,jests,bestvar,nodexts,
     &		tmissts,ncts,fill,mfixall,code,errts,dgini,tgini,
     &		tnodewt,xxts,jxxts,qxxts,nodexxts,ndpoints,nprox,jprox,npcase,
     &		kprox,noutlier,nearts,outts,wr,ncpts,touts,ndpointxxs,
     &		idataout,imax,near,iout,look,lookcls)
c
		stop
c
	endif
	END
c
c	=======================================================
c	**************  END MAIN ******************************
c	=======================================================
c
c	=======================================================
c	**************  SUBROUTINES AND FUNCTIONS  ************
c	=======================================================
c
c	-------------------------------------------------------
	subroutine runforest(mdim,ntest,nclass,maxcat,nrnodes,
     &	labelts,jbt,clts,xts,xbestsplit,qts,treemap,nbestcat,
     &	nodestatus,cat,nodeclass,jts,jests,bestvar,nodexts,
     &	tmissts,ncts,fill,mfixall,code,errts,dgini,tgini,
     &	tnodewt,xxts,jxxts,qxxts,nodexxts,ndpoints,nprox,jprox,npcase,
     &	kprox,noutlier,nearts,outts,wr,ncpts,touts,ndpointxxs,
     &	idataout,imax,near,iout,look,lookcls)
c
c	reads a forest file and runs new data through it
c
	real xts(mdim,ntest),xbestsplit(nrnodes),tmissts(nclass),
     &	qts(nclass,ntest),fill(mdim,nclass),dgini(nrnodes),
     &	tgini(mdim),tnodewt(nrnodes),xxts(mdim,ntest*nclass),
     &	qxxts(nclass,nclass*ntest),outts(ntest),wr(nclass),
     &	touts(nearts)
c
	integer treemap(2,nrnodes),nodestatus(nrnodes),
     &	cat(mdim),nodeclass(nrnodes),jests(ntest),
     &	bestvar(nrnodes),jts(ntest),clts(ntest),
     &	nodexts(ntest),ncts(nclass),jxxts(nclass*ntest),
     &	nodexxts(nclass*ntest),ndpoints(ntest,jbt),
     &	jprox(near,jbt),npcase(near,jbt),kpr(near),kprox(near,jbt),
     &	ncpts(nclass),nbestcat(maxcat,nrnodes),
     &	ndpointxxs(nclass*ntest,jbt),imax(ntest),iout(nearts)
c
	integer mdim, ntest, nclass, maxcat, nprox, nrnodes, labelts, 
     &	jbt, mfixall, code, noutlier, nearts, near, look, lookcls, 
     &	k, nts, idataout
	real errts
	integer m, j, n, jb, ndbigtree, idummy
	call zermr(qts,nclass,ntest)
	if(labelts.eq.0.and.mfixall.eq.2) then
		call zermr(qxxts,nclass,nclass*ntest) 
	endif
	read(1,*) (cat(m),m=1,mdim)
c
	if(mfixall.eq.2) then
		read(3,*) ((fill(m,j),m=1,mdim),j=1,nclass)
c		fast fix on the test data - different method if labeled
		if(labelts.eq.0) then
			call fixtestunlab(xxts,xts,mdim,nclass,ntest,fill,code)
		else
			call fixtestlab(xts,mdim,nclass,ntest,fill,code,clts)
		endif
	endif
	if(labelts.eq.1) then
		do n=1,ntest
			ncts(clts(n))=ncts(clts(n))+1
		enddo
	endif
c
c	START DOWN FOREST
c
	do jb=1,jbt
		read(1,*) ndbigtree
		do n=1,ndbigtree
			read(1,*) idummy,nodestatus(n),bestvar(n),
     &			treemap(1,n),treemap(2,n),nodeclass(n),
     &			xbestsplit(n),tnodewt(n),(nbestcat(k,n),k=1,maxcat)
    		enddo
		if(mfixall.eq.0.or.labelts.eq.1) then
			call testreebag(xts,ntest,mdim,treemap,nodestatus,
     & 			xbestsplit,bestvar,nodeclass,nrnodes,ndbigtree,
     &			cat,jts,nodexts,maxcat,nbestcat,dgini,tgini)
c
    			do n=1,ntest
				qts(jts(n),n)=qts(jts(n),n)+tnodewt(nodexts(n))
			enddo
			if(noutlier.eq.1) then
				do n=1,ntest
					ndpoints(n,jb)=nodexts(n)
				enddo
			endif
		endif !ml
		if(mfixall.eq.2.and.labelts.eq.0)then
c			unlabelled test set with missing values - use xxts
			nts=nclass*ntest
			call testreebag(xxts,nts,mdim,treemap,nodestatus,
     & 			xbestsplit,bestvar,nodeclass,nrnodes,ndbigtree,
     &			cat,jxxts,nodexxts,maxcat,nbestcat,dgini,tgini)
c
			do n=1,nclass*ntest
				k=(n/nclass)+1
				qxxts(jxxts(n),n)=qxxts(jxxts(n),n)+tnodewt(nodexxts(k))
			enddo 
			if(noutlier.eq.2) then
				do n=1,nclass*ntest
					ndpointxxs(n,jb)=nodexxts(n)
				enddo
			endif
		endif !ml
		if(labelts.eq.1) then
			if(mod(jb,look).eq.0.and.jb.le.jbt) then
				call comperrts(qts,clts,ntest,nclass,errts,
     & 				tmissts,ncts,jests)
				if(lookcls.eq.1) then
					write(*,'(i8,100f10.2)') 
     &					jb,100*errts,(100*tmissts(j), j=1,nclass)
				else
					write(*,'(i8,2f10.2)') jb,100*errts
				endif
			endif
		endif
c
	enddo !jb
c
	if(idataout.eq.2) then
		if(mfixall.eq.2.and.labelts.eq.0.and.ntest.gt.0) then
			call reducemult(wr,nclass,ntest,qxxts,jests,qts,
     &			xxts,xts,mdim,imax)
		endif
		if(labelts.eq.1) then
			do n=1,ntest
				write(7,'(3i5,1000f10.3)') n,clts(n),jests(n),
     & 				(qts(j,n), j=1,nclass),(xts(m,n),m=1,mdim)
			enddo
		else
			do n=1,ntest
				write(7,'(3i5,1000f10.3)') n,jests(n),
     & 				(qts(j,n), j=1,nclass),(xts(m,n),m=1,mdim)
			enddo
		endif
	endif
	close(7)
c
	if(noutlier.eq.2.and.labelts.eq.0.and.mfixall.eq.2) then
		do n=1,ntest
			do jb=1,jbt
				ndpoints(n,jb)=ndpointxxs((n-1)*nclass+imax(n),jb)
			enddo
		enddo
	endif
	if (noutlier.eq.2) then
		do jb=1,jbt
			do n=1,near
				read(4,*) kprox(n,jb),jprox(n,jb)
			enddo
		enddo
		call locateouts(clts,nearts,labelts,
     &		outts,jests,ncpts,ndpoints,jprox,npcase,kpr,kprox,
     &		near,jbt,touts,iout)
		do n=1,nearts
			write(15,'(2i5,f10.3)') n,jests(n),amax1(outts(n),0.0)
		enddo
		close(15)
	endif
	end
c
c	-------------------------------------------------------
	subroutine makea(x,mdim,nsample,cat,isort,v,a,b,mbax,mdimt,msm)
c
	real x(mdim,nsample),v(nsample)
	integer cat(mdim),isort(nsample),a(mdim,nsample),
     &	b(mdim,nsample),mbax(mdim),msm(mdim)
	integer mdim, nsample, mdimt
	integer k, mvar, n, n1, n2, ncat
c
c	submakea constructs the mdim x nsample integer array a.
c	If there are less than 32,000 cases, this can be declared 
c	integer*2, otherwise integer*4. For each numerical variable 
c	with values x(m,n), n=1,...,nsample, the x-values are sorted 
c	from lowest to highest. Denote these by xs(m,n). 
c	Then a(m,n) is the case  number in which 
c	xs(m,n) occurs. The b matrix is also contructed here. if the mth 
c	variable is categorical, then a(m,n) is the category of the nth 
c	case number.  
c
c	input: x, cat
c	output: a, b
c	work: v, isort
	do k=1,mdimt
	  mvar=msm(k)
	  if (cat(mvar).eq.1) then
		do n=1,nsample
			v(n)=x(mvar,n)
			isort(n)=n
		enddo
		call quicksort(v,isort,1,nsample,nsample)
c		this sorts the v(n) in ascending order. isort(n) is the 
c		case number of that v(n) nth from the lowest (assume the original
c		case numbers are 1,2,...).  
		do n=1,nsample-1
			n1=isort(n)
			n2=isort(n+1)
			a(mvar,n)=n1
			if(n.eq.1) b(mvar,n1)=1
			if (v(n).lt.v(n+1)) then
				b(mvar,n2)=b(mvar,n1)+1
			else
				b(mvar,n2)=b(mvar,n1)
			endif
		enddo
		a(mvar,nsample)=isort(nsample)
		mbax(mvar)=b(mvar,isort(nsample))
	  else
		do ncat=1,nsample
			a(mvar,ncat)=nint(x(mvar,ncat))
		enddo
	  endif
	enddo
c
	end
c
c	-------------------------------------------------------
	subroutine moda(a,nuse,nsample,mdim,cat,maxcat,
     &	ncase,jin,mdimt,msm)
c
	integer a(mdim,nsample),cat(mdim),jin(nsample),
     &	ncase(nsample),msm(mdim)
	integer nuse, nsample, mdim, mdimt, maxcat
	integer n, jj, m, k, nt, j 
	nuse=0
	do n=1,nsample
		if(jin(n).ge.1) nuse=nuse+1
	enddo
	do jj=1,mdimt
		m=msm(jj)
		k=1
		nt=1
		if(cat(m).eq.1) then
			do n=1,nsample
				if(jin(a(m,k)).ge.1) then
					a(m,nt)=a(m,k)
					k=k+1
				else
					do j=1,nsample-k
						if(jin(a(m,k+j)).ge.1) then
							a(m,nt)=a(m,k+j)
							k=k+j+1
							goto 28
						endif
					enddo
				endif
28				continue
				nt=nt+1
				if(nt.gt.nuse) goto 37
			enddo
37			continue
		endif
	enddo
	if(maxcat.gt.1) then
		k=1
		nt=1
		do n=1,nsample
			if(jin(k).ge.1) then
				ncase(nt)=k
				k=k+1
			else
				do j=1,nsample-k
					if(jin(k+j).ge.1) then
						ncase(nt)=k+j
						k=k+j+1
						goto 58
					endif
				enddo
			endif
58			continue
			nt=nt+1
			if(nt.gt.nuse) goto 85
		enddo
85		continue
	endif
	end
c
c	-------------------------------------------------------
	subroutine buildtree(a,b,cl,cat,mdim,nsample,nclass,treemap,
     &  bestvar,bestsplit,bestsplitnext,dgini, nodestatus,nodepop,
     &	nodestart,classpop,tclasspop,tclasscat,ta,nrnodes,
     &	idmove,ndsize,ncase,parent,mtry,nodeclass,ndbigtree,
     &	win,wr,wl,nuse,kcat,ncatsplit,xc,dn,cp,cm,maxcat,
     &	nbestcat,msm,mdimt)
c
c	Buildtree consists of repeated calls to findbestsplit and movedata.
c	Findbestsplit does just that--it finds the best split of the current 
c	node.  Movedata moves the data in the split node right and left so that the data
c	corresponding to each child node is contiguous.  
c
c	The buildtree bookkeeping is different from that in Friedman's original CART program. 
c	ncur is the total number of nodes to date.  nodestatus(k)=1 if the kth node has been split.
c	nodestatus(k)=2 if the node exists but has not yet been split, and =-1 of the node is
c	terminal.  A node is terminal if its size is below a threshold value, or if it is all
c	one class, or if all the x-values are equal.  If the current node k is split, then its
c	children are numbered ncur+1 (left), and ncur+2(right), ncur increases to ncur+2 and
c	the next node to be split is numbered k+1.  When no more nodes can be split, buildtree
c	returns to the main program.
c
	integer cl(nsample),cat(mdim),ncatsplit(maxcat),
     &  treemap(2,nrnodes),bestvar(nrnodes),nodeclass(nrnodes),
     &	bestsplit(nrnodes), nodestatus(nrnodes),ta(nsample),
     &  nodepop(nrnodes),nodestart(nrnodes),idmove(nsample),
     &	bestsplitnext(nrnodes),ncase(nsample),parent(nrnodes),
     &	kcat(maxcat),msm(mdim)
    	
	integer a(mdim,nsample),b(mdim,nsample)
	real tclasspop(nclass),classpop(nclass,nrnodes),
     &  tclasscat(nclass,maxcat),win(nsample),wr(nclass),
     &	wl(nclass),dgini(nrnodes),xc(maxcat),dn(maxcat),
     &	cp(maxcat),cm(maxcat)
	integer mdim, nsample, nclass, nrnodes, ndsize, mtry, ndbigtree, nuse
	integer maxcat, mdimt, j, ncur, kbuild, ndstart, ndend, jstat
	integer msplit, nbest
	integer lcat, i, kn, k, n, ndendl, nc
	real decsplit, popt1, popt2, pp
	integer nbestcat(maxcat,nrnodes)
c
	call zerv(nodestatus,nrnodes)
	call zerv(nodestart,nrnodes)
	call zerv(nodepop,nrnodes)
	call zermr(classpop,nclass,nrnodes)
	call zerm(nbestcat,maxcat,nrnodes)
	do j=1,nclass
		classpop(j,1)=tclasspop(j)
	enddo
c
	ncur=1
	nodestart(1)=1
	nodepop(1)=nuse
	nodestatus(1)=2
c	start main loop
c
	do 30 kbuild=1,nrnodes
		if (kbuild.gt.ncur) goto 50
		if (nodestatus(kbuild).ne.2) goto 30
c		initialize for next call to findbestsplit
c
		ndstart=nodestart(kbuild)
		ndend=ndstart+nodepop(kbuild)-1
		do j=1,nclass
			tclasspop(j)=classpop(j,kbuild)
		enddo
		jstat=0
		call findbestsplit(a,b,cl,mdim,nsample,nclass,cat,
     &		ndstart,ndend,tclasspop,tclasscat,msplit,decsplit,nbest,
     &		ncase,jstat,mtry,win,wr,wl,kcat,ncatsplit,xc,dn,
     &		cp,cm,maxcat,msm,mdimt)
c
		if(jstat.eq.1) then
			nodestatus(kbuild)=-1
			goto 30
		else
			bestvar(kbuild)=msplit
			dgini(kbuild)=decsplit
			if (cat(msplit).eq.1) then
				bestsplit(kbuild)=a(msplit,nbest)
				bestsplitnext(kbuild)=a(msplit,nbest+1)
			else
				lcat=cat(msplit)
				do i=1,lcat
					nbestcat(i,kbuild)=ncatsplit(i)
				enddo
			endif
		endif
		call movedata(a,ta,mdim,nsample,ndstart,ndend,idmove,ncase,
     &		msplit,cat,nbest,ndendl,ncatsplit,maxcat,mdimt,msm)
c
c		leftnode no.= ncur+1, rightnode no. = ncur+2.
c
		nodepop(ncur+1)=ndendl-ndstart+1
		nodepop(ncur+2)=ndend-ndendl
		nodestart(ncur+1)=ndstart
		nodestart(ncur+2)=ndendl+1
c
c
c		find class populations in both nodes
		do n=ndstart,ndendl
			nc=ncase(n)
			j=cl(nc)
			classpop(j,ncur+1)=classpop(j,ncur+1)+win(nc)
		enddo
		do n=ndendl+1,ndend
			nc=ncase(n)
			j=cl(nc)
			classpop(j,ncur+2)=classpop(j,ncur+2)+win(nc)
		enddo
c
c		check on nodestatus
c
		nodestatus(ncur+1)=2
		nodestatus(ncur+2)=2
		if (nodepop(ncur+1).le.ndsize) nodestatus(ncur+1)=-1
		if (nodepop(ncur+2).le.ndsize) nodestatus(ncur+2)=-1
		popt1=0
		popt2=0
		do j=1,nclass
			popt1=popt1+classpop(j,ncur+1)
			popt2=popt2+classpop(j,ncur+2)
		enddo
		do j=1,nclass
			if (abs(classpop(j,ncur+1)-popt1).lt.8.232D-11) nodestatus(ncur+1)=-1
			if (abs(classpop(j,ncur+2)-popt2).lt.8.232D-11) nodestatus(ncur+2)=-1
		enddo
		treemap(1,kbuild)=ncur+1
		treemap(2,kbuild)=ncur+2
		parent(ncur+1)=kbuild
		parent(ncur+2)=kbuild
		nodestatus(kbuild)=1
		ncur=ncur+2
		if (ncur.ge.nrnodes) goto 50
30	continue
50	continue
	ndbigtree=nrnodes
	do k=nrnodes,1,-1
		if (nodestatus(k).eq.0) ndbigtree=ndbigtree-1
		if (nodestatus(k).eq.2) nodestatus(k)=-1
	enddo
	do kn=1,ndbigtree
		if(nodestatus(kn).eq.-1) then
			pp=0
			do j=1,nclass
				if(classpop(j,kn).gt.pp) then
					nodeclass(kn)=j
					pp=classpop(j,kn)
				endif
			enddo
		endif
	enddo
	end
c
c	-------------------------------------------------------
	subroutine findbestsplit(a,b,cl,mdim,nsample,nclass,cat,
     &  ndstart,ndend,tclasspop,tclasscat,msplit,decsplit,nbest,
     &  ncase,jstat,mtry,win,wr,wl,kcat,ncatsplit,xc,dn,
     &	cp,cm,maxcat,msm,mdimt)
c
c	For the best split, msplit is the variable split on. decsplit is the dec. in impurity. 
c	If msplit is numerical, nsplit is the case number of value of msplit split on,
c	and nsplitnext is the case number of the next larger value of msplit.  If msplit is
c	categorical, then nsplit is the coding into an integer of the categories going left.
c
	integer a(mdim,nsample),b(mdim,nsample)
c
	integer cl(nsample),cat(mdim), ncase(nsample),msm(mdimt),
     &	kcat(maxcat),ncatsplit(maxcat),icat(32)
c
	integer mdim, nsample, nclass, ndstart, ndend, mdimt, mtry,
     &	msplit, nbest, jstat, maxcat, i, j, k, l, mv, mvar, nc, 
     &	lcat, nnz, nhit, ncatsp, nsp
c
	real tclasspop(nclass),tclasscat(nclass,maxcat),
     &	win(nsample), wr(nclass),wl(nclass),xc(maxcat),
     &	dn(maxcat),cp(maxcat),cm(maxcat)
c
	real decsplit, pno, pdo, rrn, rrd, rln, rld, u,
     &	crit0, critmax, crit, su
c
	real rrand
c
	external unpack
c
c	compute initial values of numerator and denominator of Gini
	pno=0
	pdo=0
	do j=1,nclass
		pno=pno+tclasspop(j)*tclasspop(j)
		pdo=pdo+tclasspop(j)
	enddo
	crit0=pno/pdo
	jstat=0
c	start main loop through variables to find best split
	critmax=-1.0e20
	do 20 mv=1,mtry
200		continue
		k=int(mdimt*rrand(1))+1
		mvar=msm(k)
		if(cat(mvar).eq.1) then
			rrn=pno
			rrd=pdo
			rln=0
			rld=0
			call zervr(wl,nclass)
			do j=1,nclass
				wr(j)=tclasspop(j)
			enddo
			do nsp=ndstart,ndend-1
				nc=a(mvar,nsp)
				u=win(nc)
				k=cl(nc)
				rln=rln+u*(2*wl(k)+u)
				rrn=rrn+u*(-2*wr(k)+u)
				rld=rld+u
				rrd=rrd-u
				wl(k)=wl(k)+u
				wr(k)=wr(k)-u
  				if (b(mvar,nc).lt.b(mvar,a(mvar,nsp+1))) then
					if(amin1(rrd,rld).gt.1.0e-5) then
						crit=(rln/rld)+(rrn/rrd)
						if (crit.gt.critmax) then
							nbest=nsp
							critmax=crit
							msplit=mvar
						endif
					endif
				endif
			enddo
 		else
c			compute the decrease in impurity given by categorical splits
			lcat=cat(mvar)
			call zermr(tclasscat,nclass,maxcat)
			do nsp=ndstart,ndend
				nc=ncase(nsp)
				l=a(mvar,ncase(nsp))
				tclasscat(cl(nc),l)=tclasscat(cl(nc),l)+win(nc)
			enddo
			nnz=0
			do i=1,lcat
				su=0
				do j=1,nclass
					su=su+tclasscat(j,i)
				enddo
				dn(i)=su
				nnz=nnz+1
			enddo
			if (nnz.eq.1) then
				critmax=-1.0e25
				goto 20
			endif
			if(lcat.gt.10.and.nclass.eq.2 )then
				call catmaxb(tclasscat,tclasspop,xc,cp,cm,kcat,nclass,
     &				lcat,maxcat,ncatsplit,critmax,pdo,nhit,dn)
				if(nhit.eq.1)then
					msplit=mvar
				endif
			endif
			if(lcat.le.10) then
				call catmax(pdo,tclasscat,tclasspop,nclass,lcat,
     & 				ncatsp,critmax,nhit,maxcat)
				if(nhit.eq.1) then
					msplit=mvar
					call unpack(lcat,ncatsp,icat)
					call zerv(ncatsplit,maxcat)
					do k=1,lcat
						ncatsplit(k)=icat(k)
					enddo
				endif
			endif
			if(lcat.gt.10.and.nclass.gt.2) then
				call catmaxap(tclasscat,kcat,nclass,lcat,
     &				maxcat,critmax,dn,nhit)
				if(nhit.eq.1) then
					msplit=mvar
					call zerv(ncatsplit,maxcat)
					do l=1,lcat
						if(kcat(l).eq.-1) ncatsplit(l)=1
					enddo
				endif
			endif
		endif !cat
20	continue
25	continue
	decsplit=critmax-crit0
	if (critmax.lt.-1.0e10) jstat=1
	end
c
c	-------------------------------------------------------
	subroutine catmaxb(tclasscat,tclasspop,xc,cp,cm,kcat,
     &	nclass,lcat,maxcat,ncatsplit,critmax,pdo,nhit,dn)
c
	real tclasscat(nclass,maxcat),tclasspop(nclass),
     &	xc(maxcat),cp(maxcat),cm(maxcat),dn(maxcat)
	real critmax, pdo
c
	integer kcat(maxcat),ncatsplit(maxcat)
	integer nclass, lcat, maxcat, nhit, l, nk, j, n, i, k
	real bestsplit, rrd, rld, rln, rrn, crit
	nhit=0
 	do l=1,lcat
		if(dn(l).gt.0) then
			xc(l)=tclasscat(1,l)/dn(l)
		else
			xc(l)=0
		endif
	enddo
	do nk=1,lcat
		kcat(nk)=nk
	enddo
	call quicksort(xc,kcat,1,lcat,maxcat)
	do j=1,nclass
		cp(j)=0
		cm(j)=tclasspop(j)
	enddo
	rrd=pdo
	rld=0
	do n=1,lcat-1
		rld=rld+dn(kcat(n))
		rrd=rrd-dn(kcat(n))
		do k=1,nclass
			cp(k)=cp(k)+tclasscat(k,kcat(n))
			cm(k)=cm(k)-tclasscat(k,kcat(n))
		enddo
		rln=0
		rrn=0
		do k=1,nclass
			rln=rln+cp(k)**2
			rrn=rrn+cm(k)**2
		enddo
		if (xc(n).lt.xc(n+1)) then
			if(amin1(rrd,rld).gt.1) then
				crit=(rln/rld)+(rrn/rrd)
				if (crit.gt.critmax) then
					critmax=crit
					bestsplit=.5*(xc(n)+xc(n+1))
					nhit=1
				endif
			endif
		endif
	enddo !n
	if (nhit.eq.1) then
		call zerv(ncatsplit,maxcat)
		do l=1,lcat
			if(dn(l).gt.0) then
				xc(l)=tclasscat(1,l)/dn(l)
			else
				xc(l)=0
			endif
		enddo
		do i=1,lcat
			if(xc(i).lt.bestsplit) ncatsplit(i)=1
		enddo
	endif
	end
c
c	-------------------------------------------------------
	subroutine catmax(pdo,tclasscat,tclasspop,nclass,lcat,
     &  ncatsp,critmax,nhit,maxcat)
c
c	this finds the best categorical split of a categorical variable
c	with lcat categories, nclass classes and tclasscat(j,l) is the number of cases in
c	class j with category value l. The method used is an exhaustive search over all
c	partitions of the category values.  For the two class problem, there is a
c	faster exact algorithm we will add later.  If lcat.ge.10, the exhaustive search
c	gets slow and there is a faster iterative algorithm we can add later.
	real tclasscat(nclass,maxcat),tclasspop(nclass),tmpclass(100)
	integer icat(32),n, lcat, nhit, l, j, ncatsp, nclass, maxcat
	real critmax, pdo, pln, pld, prn, tdec
c
	external unpack
	nhit=0
	do n=1,(2**(lcat-1))-1
		call unpack(lcat,n,icat)
		do j=1,nclass
			tmpclass(j)=0
			do l=1,lcat
				if(icat(l).eq.1) then
					tmpclass(j)=tmpclass(j)+tclasscat(j,l)
				endif
			enddo
		enddo
		pln=0
		pld=0
		do j=1,nclass
			pln=pln+tmpclass(j)*tmpclass(j)
			pld=pld+tmpclass(j)
		enddo
		prn=0
		do j=1,nclass
			tmpclass(j)=tclasspop(j)-tmpclass(j)
			prn=prn+tmpclass(j)*tmpclass(j)
		enddo
		tdec=(pln/pld)+(prn/(pdo-pld))
		if (tdec.gt.critmax) then
			critmax=tdec
			ncatsp=n
			nhit=1
		endif
	enddo
	end
c
c	-------------------------------------------------------
	subroutine catmaxap(tclasscat,kcat,
     &	nclass,lcat,maxcat,critmax,dn,nhit)
c
c	this returns those categories going left in the best split. 
c	This is coded into a vector of length maxcat (see under catmax below for details). 
c
	real tclasscat(nclass,maxcat),dn(maxcat),ncr0(100),
     &	ncl(100),ncr(100),ncl0(100),ntr,ntl,ntl0,ntr0
c
	integer kcat(maxcat)
c
	real critmax, crit, critnew
	integer nclass, lcat, maxcat, nhit, l, j, k, nchange
c
	real rrand
	nhit=0
	do l=1,lcat
		if(rrand(1).le..5) then
			kcat(l)=-1
		else
			kcat(l)=1
		endif
	enddo
	do j=1,nclass
		ncl(j)=0
		ncr(j)=0
		ntl=0
		ntr=0
	enddo
	do l=1,lcat
		do j=1,nclass
			if(kcat(l).eq.-1) then
				ncl(j)=ncl(j)+tclasscat(j,l)
				ntl=ntl+dn(l)
			else
				ncr(j)=ncr(j)+tclasscat(j,l)
				ntr=ntr+dn(l)
			endif
		enddo
	enddo
	crit=0
	do j=1,nclass
		crit=crit+(ncl(j)*ncl(j)/ntl)+(ncr(j)*ncr(j)/ntr)
	enddo
	do k=1,1000
		nchange=0
		do l=1,lcat
			do j=1,nclass
				ncl0(j)=ncl(j)
				ncr0(j)=ncr(j)
				ntl0=ntl
				ntr0=ntr
				ncl(j)=ncl(j)+kcat(l)*tclasscat(j,l)
				ncr(j)=ncr(j)-kcat(l)*tclasscat(j,l)
				ntl=ntl+kcat(l)*dn(l)
				ntr=ntr-kcat(l)*dn(l)
			enddo
			critnew=0
			do j=1,nclass
				critnew=critnew+(ncl(j)*ncl(j)/ntl)+(ncr(j)*ncr(j)/ntr)
			enddo
			if(critnew.gt.critmax) then
				critmax=critnew
				nchange=nchange+1
				kcat(l)=-kcat(l)
				nhit=1
			else
				do j=1,nclass
					ncl(j)=ncl0(j)
					ncr(j)=ncr0(j)
					ntl=ntl0
					ntr=ntr0
				enddo
			endif
		enddo !l
		if(nchange.eq.0) goto 101
	enddo !k
101	continue
	end
c
c	-------------------------------------------------------
	subroutine movedata(a,ta,mdim,nsample,ndstart,ndend,idmove,
     &  ncase,msplit,cat,nbest,ndendl,ncatsplit,maxcat,mdimt,msm)
c
c	movedata is the heart of the buildtree construction. 
c	Based on the best split the data corresponding to the 
c	current node is moved to the left if it belongs to the 
c	left child and right if it belongs to the right child.
	integer a(mdim,nsample),ta(nsample),idmove(nsample),
     &  ncase(ndend),cat(mdim),ncatsplit(maxcat),msm(mdimt)
	integer mdim, ndstart, ndend, nsample, msplit, nbest, ndendl
	integer maxcat, mdimt, nsp, nc, ms, msh, k, n, ih
c	compute idmove=indicator of case nos. going left
c
	if(cat(msplit).eq.1) then
c		continuous variable
		do nsp=ndstart,nbest
			nc=a(msplit,nsp)
			idmove(nc)=1
		enddo
		do nsp=nbest+1,ndend
			nc=a(msplit,nsp)
			idmove(nc)=0
		enddo
		ndendl=nbest
	else
c		categorical variable
		ndendl=ndstart-1
		do nsp=ndstart,ndend
			nc=ncase(nsp)
			if(ncatsplit(a(msplit,nc)).eq.1) then
				idmove(nc)=1
				ndendl=ndendl+1
			else
				idmove(nc)=0
			endif
		enddo
	endif
c	shift case nos. right and left 
	do ms=1,mdimt
		msh=msm(ms)
		k=ndstart-1
		if(cat(msh).eq.1) then
c			continuous variable
			do n=ndstart,ndend
				ih=a(msh,n)
				if(idmove(ih).eq.1) then
					k=k+1
					ta(k)=a(msh,n)
				endif
			enddo
			do n=ndstart,ndend
				ih=a(msh,n)
				if(idmove(ih).eq.0) then 
					k=k+1
					ta(k)=a(msh,n)
				endif
			enddo
		else
c			categorical variable
			do n=ndstart,ndend
				ih=ncase(n)
				if(idmove(ih).eq.1) then
					k=k+1
					ta(k)=a(msh,ih)
				endif
			enddo
			do n=ndstart,ndend
				ih=ncase(n)
				if(idmove(ih).eq.0) then 
					k=k+1
					ta(k)=a(msh,ih)
				endif
			enddo
		endif
		do k=ndstart,ndend
			a(msh,k)=ta(k)
		enddo
	enddo
c
c	compute case nos. for right and left nodes.
	if(cat(msplit).eq.1) then
		do n=ndstart,ndend
			ncase(n)=a(msplit,n)
		enddo
	else
		k=ndstart-1
		do n=ndstart,ndend
			if(idmove(ncase(n)).eq.1) then
				k=k+1
				ta(k)=ncase(n)
			endif
		enddo
		do n=ndstart,ndend
			if(idmove(ncase(n)).eq.0) then
				k=k+1
				ta(k)=ncase(n)
			endif
		enddo
		do k=ndstart,ndend
			ncase(k)=ta(k)
		enddo
	endif
	end
c
c	-------------------------------------------------------
	subroutine xtranslate(x,mdim,nrnodes,nsample,bestvar,
     &  bestsplit,bestsplitnext,xbestsplit,nodestatus,cat,
     &	ndbigtree)
c
c	xtranslate takes the splits on numerical variables and translates them
c	back into x-values.  It also unpacks each categorical split into a 32-dimensional 
c	vector with components of zero or one--a one indicates that the corresponding
c	category goes left in the split.
c
	integer cat(mdim),bestvar(nrnodes),bestsplitnext(nrnodes),
     &	nodestatus(nrnodes),bestsplit(nrnodes)
	real x(mdim,nsample),xbestsplit(nrnodes)
	integer mdim, nrnodes, nsample, ndbigtree, k, m
c
	do k=1,ndbigtree
		if (nodestatus(k).eq.1) then
			m=bestvar(k)
			if (cat(m).eq.1) then
				xbestsplit(k)=(x(m,bestsplit(k))+
     &				x(m,bestsplitnext(k)))/2
			else
				xbestsplit(k)=real(bestsplit(k))
			endif
		endif
	enddo
	end
c
c	-------------------------------------------------------
	subroutine getweights(x,nsample,mdim,treemap,nodestatus,
     &  xbestsplit,bestvar,nrnodes,ndbigtree,
     &	cat,maxcat,nbestcat,jin,win,tw,tn,tnodewt)
c
	real x(mdim,nsample),xbestsplit(nrnodes),
     &	win(nsample),tw(ndbigtree),tn(ndbigtree),tnodewt(ndbigtree)
	integer nsample, mdim, nrnodes, ndbigtree, maxcat
	integer treemap(2,nrnodes),bestvar(nrnodes),
     &  cat(mdim),nodestatus(nrnodes),jin(nsample)
c
	integer nbestcat(maxcat,nrnodes)
	integer jcat, n, kt, k, m
	call zervr(tw,ndbigtree)
	call zervr(tn,ndbigtree)
	do n=1,nsample
		if(jin(n).ge.1) then
			kt=1
			do k=1,ndbigtree
				if (nodestatus(kt).eq.-1) then
					tw(kt)=tw(kt)+win(n)
					tn(kt)=tn(kt)+jin(n)
					goto 100
				endif
				m=bestvar(kt)
				if (cat(m).eq.1) then
					if (x(m,n).le.xbestsplit(kt)) then 
						kt=treemap(1,kt)
					else
						kt=treemap(2,kt)
					endif
				else
					jcat=nint(x(m,n))
					if (nbestcat(jcat,kt).eq.1) then
						kt=treemap(1,kt)
					else
						kt=treemap(2,kt)
					endif
				endif
			enddo
100			continue
		endif
	enddo
	do n=1,ndbigtree
		if(nodestatus(n).eq.-1) tnodewt(n)=tw(n)/tn(n)
	enddo
	end
c
c	-------------------------------------------------------
	subroutine testreebag(x,nts,mdim,treemap,nodestatus,
     &  xbestsplit,bestvar,nodeclass,nrnodes,ndbigtree,
     &	cat,jtr,nodextr,maxcat,nbestcat,dgini,tgini)
c
c	predicts the class of all objects in x
c
c	input:
	real x(mdim,nts),xbestsplit(nrnodes),dgini(nrnodes)
	integer treemap(2,nrnodes),bestvar(nrnodes),
     &  nodeclass(nrnodes),cat(mdim),nodestatus(nrnodes),
     &	nbestcat(maxcat,nrnodes)
	integer nts, mdim, nrnodes, ndbigtree, maxcat
c
c	output:
	real tgini(mdim)
	integer nodextr(nts),jtr(nts)
c
c	local
	integer n, kt, k, m, jcat
c
	call zerv(jtr,nts)
	call zerv(nodextr,nts)
	call zervr(tgini,mdim)
	n=1
903	kt=1
	do k=1,ndbigtree
		if (nodestatus(kt).eq.-1) then
			jtr(n)=nodeclass(kt)
			nodextr(n)=kt
			goto 100
		endif
		m=bestvar(kt)
		tgini(m)=tgini(m)+dgini(kt)
		if (cat(m).eq.1) then
			if (x(m,n).le.xbestsplit(kt)) then 
				kt=treemap(1,kt)
			else
				kt=treemap(2,kt)
			endif
		endif
		if(cat(m).gt.1) then
			jcat=nint(x(m,n))
			if(jcat.lt.1) write(*,*) m,n,jcat
			if (nbestcat(jcat,kt).eq.1) then
				kt=treemap(1,kt)
			else
				kt=treemap(2,kt)
			endif
		endif
	enddo !k
100	n=n+1
	if(n.le.nts) goto 903
	do m=1,mdim
		tgini(m)=tgini(m)/nts
	enddo
	end
c
c	-------------------------------------------------------
	subroutine testreeimp(x,nsample,mdim,joob,pjoob,nout,mr,
     &	treemap,nodestatus,xbestsplit,bestvar,nodeclass,nrnodes,
     &	ndbigtree,cat,jvr,nodexvr,maxcat,nbestcat)
c
c	predicts the class of out-of-bag-cases for variable importance
c	also computes nodex
c
c	input:
	real x(mdim,nsample),xbestsplit(nrnodes)
	integer joob(nout),pjoob(nout),
     &	treemap(2,nrnodes),bestvar(nrnodes),
     &  nodeclass(nrnodes),cat(mdim),nodestatus(nrnodes),
     &	nbestcat(maxcat,nrnodes)
	integer nsample, mdim, nout, mr,
     &	nrnodes, ndbigtree, maxcat
c
c	output:
	integer nodexvr(nout),jvr(nout)
c
c	local
	integer n, kt, k, m, jcat
	real xmn
c
	call permobmr(joob,pjoob,nout)
	n=1
904	kt=1
	do k=1,ndbigtree
		if (nodestatus(kt).eq.-1) then
			jvr(n)=nodeclass(kt)
			nodexvr(n)=kt
			goto 100
		endif
		m=bestvar(kt)
		if(m.eq.mr) then
			xmn = x(m,pjoob(n)) ! permuted value
		else
			xmn = x(m,joob(n))
		endif
		if (cat(m).eq.1) then
			if (xmn.le.xbestsplit(kt)) then 
				kt=treemap(1,kt)
			else
				kt=treemap(2,kt)
			endif
		endif
		if(cat(m).gt.1) then
			jcat=nint(xmn)
			if(jcat.lt.1) write(*,*) m,n,jcat
			if (nbestcat(jcat,kt).eq.1) then
				kt=treemap(1,kt)
			else
				kt=treemap(2,kt)
			endif
		endif
	enddo !k
100	n=n+1
	if(n.le.nout) goto 904
	end
c
c	-------------------------------------------------------
	subroutine permobmr(joob,pjoob,nout)
c
c	randomly permute the elements of joob and put them in pjoob
c
c	input:
	integer joob(nout),nout
c	output:
	integer pjoob(nout)
c	local:
	integer j, k, jt
	real rnd,rrand
c
	do j=1,nout
		pjoob(j) = joob(j)
	enddo 
	j=nout
11	rnd = rrand(1)
	k=int(j*rnd)+1
c	switch j and k
	jt=pjoob(j)
	pjoob(j)=pjoob(k)
	pjoob(k)=jt
	j=j-1
	if(j.gt.1) go to 11
	end
c
c	-------------------------------------------------------
	subroutine comperrtr(q,cl,nsample,nclass,errtr,
     &  tmiss,nc,jest,out)
c
	integer cl(nsample),nc(nclass),jest(nsample),out(nsample)
	real q(nclass,nsample),tmiss(nclass), errtr
	integer nsample, nclass
	real cmax, ctemp
	integer n, j, jmax
	call zervr(tmiss,nclass)
	errtr=0
	do n=1,nsample
		cmax=0
		do j=1,nclass
			if(out(n).gt.0) ctemp=q(j,n)/out(n)
			if(ctemp.gt.cmax) then
				jmax=j
				cmax=ctemp
			endif
		enddo
		jest(n)=jmax
		if (jmax.ne.cl(n)) then
			tmiss(cl(n))=tmiss(cl(n))+1
			errtr=errtr+1
		endif
	enddo
	errtr=errtr/nsample
	do j=1,nclass
		tmiss(j)=tmiss(j)/nc(j)
	enddo
	end
c
c	-------------------------------------------------------
	subroutine comperrts(qts,clts,ntest,nclass,errts,
     &  tmissts,ncts,jests)
c
c	note: only called if there is a labeled test set
c
	integer clts(ntest),ncts(nclass),jests(ntest)
	real qts(nclass,ntest),tmissts(nclass)
	integer ntest, nclass
	real errts, cmax
	integer n, j, jmax
	call zervr(tmissts,nclass)
	errts=0
	do n=1,ntest
		cmax=0
		do j=1,nclass
			if (qts(j,n).gt.cmax) then
				jmax=j
				cmax=qts(j,n)
			endif
		enddo
		jests(n)=jmax
		if (jmax.ne.clts(n)) then
			tmissts(clts(n))=tmissts(clts(n))+1
			errts=errts+1
		endif
	enddo
	errts=errts/ntest
	do j=1,nclass
		tmissts(j)=tmissts(j)/ncts(j)
	enddo
	end
c
c	-------------------------------------------------------
	subroutine createclass(x,cl,ns,nsample,mdim)
c
	real x(mdim,nsample)
	integer cl(nsample)
	integer ns, nsample, mdim
	real rrand
	integer n, m, k
c
	do n=1,ns
		cl(n)=1
	enddo
	do n=ns+1,nsample
		cl(n)=2
	enddo
	do n=ns+1,nsample
		do m=1,mdim
			k=int(rrand(1)*ns)+1
			x(m,n)=x(m,k)
		enddo
	enddo
	end
c
c	-------------------------------------------------------
	subroutine finishimp(mdim,sqsd,avimp,signif,zscore,
     &	jbt,mdimt,msm)
c
	real sqsd(mdim),avimp(mdim),zscore(mdim),signif(mdim)
   	integer msm(mdim)
	integer mdim, mdimt, k, m1, jbt
	real v, av, se
	real erfcc
	do k=1,mdimt
		m1=msm(k)
		avimp(m1)=avimp(m1)/jbt
		av=avimp(m1)
		se=(sqsd(m1)/jbt)-av*av
		se=sqrt(se/jbt)
		if(se.gt.0.0) then
			zscore(m1)=avimp(m1)/se
			v=zscore(m1)
			signif(m1)=erfcc(v)
		else
			zscore(m1)=-5
			signif(m1)=1
		endif
	enddo
	end
c
c	-------------------------------------------------------
	subroutine compinteract(votes,effect,msm,mdim,mdimt,
     &	jbt,g,iv,irnk,hist,teffect)
c
	real votes(mdim,jbt),effect(mdim,mdim),g(mdim),
     &	hist(0:mdim,mdim),teffect(mdim,mdim)
	integer msm(mdimt),iv(mdim),irnk(mdim,jbt)
	integer mdim, mdimt, jbt, jb, i, nt, irk, j, ii, jj, ij, mmin, m, k
	real gmin, rcor
	do jb=1,jbt
		nt=0
		call zerv(iv,mdim)
		do i=1,mdimt
			m=msm(i)
			g(m)=votes(m,jb)
			if(abs(g(m)).lt.8.232D-11) then
				irnk(m,jb)=0
				iv(m)=1
				nt=nt+1
			endif
		enddo
		irk=0
		do j=1,8000
			gmin=10000
			do i=1,mdimt
				m=msm(i)
				if(iv(m).eq.0.and.g(m).lt.gmin) then
					gmin=g(m)
					mmin=m
				endif
			enddo
			iv(mmin)=1
			irk=irk+1
			irnk(mmin,jb)=irk
			nt=nt+1
			if(nt.ge.mdimt) goto 79
		enddo !j
79		continue
	enddo !jb
	do j=0,mdimt
		do i=1,mdimt
			m=msm(i)
			hist(j,m)=0
		enddo
	enddo
	do i=1,mdimt
		m=msm(i)
		do jb=1,jbt
			hist(irnk(m,jb),m)=hist(irnk(m,jb),m)+1
		enddo
		do j=0,mdimt
			hist(j,m)=hist(j,m)/jbt
		enddo
	enddo !m
	do j=0,mdimt
		!write(*,'(14f8.2)') (hist(j,msm(i)),i=1,10)
	enddo
	call zermr(effect,mdim,mdim)
	do i=1,mdimt
		do j=1,mdimt
			m=msm(i)
			k=msm(j)
			do jb=1,jbt
				effect(m,k)=effect(m,k)+iabs(irnk(m,jb)-irnk(k,jb))
			enddo
			effect(m,k)=effect(m,k)/jbt
		enddo
	enddo
	call zermr(teffect,mdim,mdim)
	do i=1,mdimt
		do j=1,mdimt
			m=msm(i)
			k=msm(j)
			do ii=0,mdimt
				do jj=0,mdimt
					teffect(m,k)=teffect(m,k)+abs(ii-jj)*hist(jj,m)*hist(ii,k)
				enddo
			enddo
			rcor=0
			do ij=1,mdimt
				rcor=rcor+hist(ij,m)*hist(ij,k)
			enddo
			teffect(m,k)=teffect(m,k)/(1-rcor)
		enddo
	enddo
	do i=1,mdimt
		do j=1,mdimt
			m=msm(i)
			k=msm(j)
			effect(m,k)=100*(effect(m,k)-teffect(m,k))
		enddo
	enddo
	end
c
c	-------------------------------------------------------
	subroutine compproto(prox,jest,fp,loz,nsample,
     &	nprot,tk,nrnn,nclass,kprot,kend)
c
	real tk(0:nprot,nsample),fp(nprot,nclass)
    	integer njt(100),jest(nsample),loz(nsample,nrnn),
     &	kprot(nprot,nclass),kend(nsample)
	double precision prox(nsample,nrnn)
	integer nsample, nprot, nclass, nrnn
	integer n, jp, ip, kp, j, lmax, i, kmax, ind, k, l
	real pmax, pmax0, tmax, aa
c
	call zermr(fp,nprot,nclass)
	call zerv(njt,100)
	do n=1,nsample
		njt(jest(n))=njt(jest(n))+1
	enddo
	do jp=1,nclass
		do n=1,nsample
			tk(0,n)=0
		enddo
		do ip=1,nprot
			do kp=1,nsample
				if(jest(kp).eq.jp) then
					pmax0=0
					pmax=0
					do n=1,nsample
						if(jest(n).eq.jp) then
							do j=1,min0(nrnn,kend(n))
								if(loz(n,j).eq.kp) then
									if(real(prox(n,j)).gt.tk(ip-1,n)) 
     &									pmax=pmax+real(prox(n,j))
								endif
							enddo
						endif
					enddo
					if(pmax.gt.pmax0) then
						kmax=kp
						pmax0=pmax
					endif
				endif
			enddo
			kprot(ip,jp)=kmax
			do n=1,nsample
				if(jest(n).eq.jp) then
					ind=0
					do i=1,min0(nrnn,kend(n))
						if(loz(n,i).eq.kprot(ip,jp)) then
							ind=1
							aa=real(prox(n,i))
							tk(ip,n)=amax1(tk(ip-1,n),aa)
						endif
					enddo
					if(ind.eq.0) tk(ip,n)=tk(ip-1,n)
				endif
			enddo !n
		enddo !ip
		do n=1,nsample
			if(jest(n).eq.jp) then
				lmax=1
				tmax=tk(1,n)
				do l=2,nprot
					if(tk(l,n).gt.tmax) then
						lmax=l
						tmax=tk(l,n)
					endif
				enddo
				fp(lmax,jp)=fp(lmax,jp)+1
			endif
		enddo
	enddo !jp
	do jp=1,nclass
		do k=1,nprot
			fp(k,jp)=100.*fp(k,jp)/njt(jp)
		enddo
	enddo
	end
c		
c	-------------------------------------------------------
	subroutine preprox(nodex,near,ncass,jpr,kpx)
	integer nodex(near),ncass(near),jpr(near),kpx(near)
	integer near, n, j, k
c
	call zerv(jpr,near)
	call zerv(kpx,near)
	do n=1,near
		ncass(n)=n
	enddo
	call quicksortint(nodex,ncass,1,near,near)
	k=1
	jpr(1)=1
	kpx(nodex(1))=1
	do j=1,near-1
		if(nodex(j+1).gt.nodex(j)) then
			k=k+1
			jpr(k)=j+1
			kpx(nodex(j+1))=k
		endif
	enddo
	jpr(k+1)=near+1
	end
c
c	-------------------------------------------------------
	subroutine comprox(kprox,ndpoint,jprox,npcase,
     &	kpr,lcase,loz,prox,near,jbt,nrnn,kend,key,inuse)
c
	double precision prox(near,nrnn)
	integer ndpoint(near,jbt),jprox(near,jbt),kpr(near),
     &	kprox(near,jbt),npcase(near,jbt),loz(near,nrnn),
     &	lcase(0:jbt),kend(near),key(near),inuse(near)
	integer near, jbt, nrnn
	integer n, j, jb, kk, jj, nn, nt, ntt, k, ntot, ncrit
c
	do n=1,near
		call zerv(kpr,near)
		call zerv(inuse,near)
		do j=0,jbt
			lcase(j)=0
		enddo
		ntt=0
		do jb=1,jbt
			k=kprox(ndpoint(n,jb),jb)
			do j=jprox(k,jb),jprox(k+1,jb)-1
				kk=npcase(j,jb)
				kpr(kk)=kpr(kk)+1
				if(inuse(kk).eq.0) then
					ntt=ntt+1
					key(ntt)=kk
					inuse(kk)=1
				endif
			enddo
		enddo
		do kk=1,ntt
			jj=kpr(key(kk))
			lcase(jj)=lcase(jj)+1
		enddo
		ntot=0
		do nn=jbt,1,-1
			ntot=ntot+lcase(nn)
			if(ntot.ge.nrnn) then
				ncrit=nn
				go to 14
			endif
		enddo
14		continue
c
		nt=0
		do k=1,ntt
			kk=key(k)
			if(kpr(kk).ge.ncrit) then
				nt=nt+1
				loz(n,nt)=kk
				prox(n,nt)=dble(kpr(kk))/jbt
				if(nt.ge.nrnn)then
					goto 24
				endif
			endif
		enddo
24		continue
		kend(n)=nrnn
	enddo !n
	end
c
c	-------------------------------------------------------
	subroutine comproxts(kprox,jprox,npcase,kpr,lcase,
     &	near,jbt,nrnn,nearts,proxts,lozts,ndpoints,kendts,
     &	key,inuse)
c
	double precision proxts(nearts,nrnn)
	integer kprox(near,jbt),jprox(near,jbt),npcase(near,jbt),kpr(near),
     &	lcase(0:jbt),inuse(near),kendts(nearts),key(near),
     &	ndpoints(nearts,jbt),lozts(nearts,nearts)
c
	integer near, nearts, jbt, nrnn
	integer n, j, jb, kk, jj, nn, nt, ntt, k, ntot, ncrit
c
	do n=1,nearts
		call zerv(kpr,near)
		call zerv(inuse,near)
		do j=0,jbt
			lcase(j)=0
		enddo
		ntt=0
		do jb=1,jbt
			k=kprox(ndpoints(n,jb),jb)
			do j=jprox(k,jb),jprox(k+1,jb)-1
				kk=npcase(j,jb)
				kpr(kk)=kpr(kk)+1
				if(inuse(kk).eq.0) then
					ntt=ntt+1
					key(ntt)=kk
					inuse(kk)=1
				endif
			enddo
		enddo
		do kk=1,ntt
			jj=kpr(key(kk))
			lcase(jj)=lcase(jj)+1
		enddo
		ntot=0
		do nn=jbt,1,-1
			ntot=ntot+lcase(nn)
			if(ntot.ge.nrnn) then
				ncrit=nn
				go to 14
			endif
		enddo
14		continue
c
		nt=0
		do k=1,ntt
			kk=key(k)
			if(kpr(kk).ge.ncrit) then
				nt=nt+1
				lozts(n,nt)=kk
				proxts(n,nt)=dble(kpr(kk))/jbt
				if(nt.ge.nrnn)then
					goto 24
				endif
			endif
		enddo
24		continue
		kendts(n)=nrnn
	enddo !n
	end
c
c	-------------------------------------------------------
	subroutine locateout(cl,near,nclass,
     &	labeltr,outtr,tout,isort,devout,rmedout,ncp,
     &	ndpoint,jprox,npcase,kpr,kprox,jbt)
c
	real outtr(near),tout(near),devout(nclass),rmedout(nclass)
	integer cl(near),isort(near), ndpoint(near,jbt),
     &	jprox(near,jbt),npcase(near,jbt),kpr(near),kprox(near,jbt),
     &	ncp(near)
c
	integer nclass, labeltr, near, jbt, n, jb, k, j, jp,
     &	jpclass, nt, jc
	real rmed, dev, rsq
c
	if(labeltr.eq.1) jpclass=nclass
	if(labeltr.eq.0) jpclass=1
	do n=1,near
		jc=min0(jpclass,cl(n))
		call zerv(kpr,near)
		do jb=1,jbt
			k=kprox(ndpoint(n,jb),jb)
			do j=jprox(k,jb),jprox(k+1,jb)-1
				if(cl(npcase(j,jb)).eq.jc) then
					kpr(npcase(j,jb))=kpr(npcase(j,jb))+1
				endif
			enddo
		enddo
		rsq=0
		do k=1,near
			if(kpr(k).gt.0) rsq=rsq+kpr(k)*kpr(k)
		enddo
		outtr(n)=near/rsq
	enddo !n
	do n=1,near
		!write(*,*) n,outtr(n)
	enddo
	do jp=1,jpclass
		nt=0
		do n=1,near
			if(cl(n).eq.jp) then
				nt=nt+1
				tout(nt)=outtr(n)
				ncp(nt)=n
			endif
		enddo
		do n=1,nt
			isort(n) = n
		enddo
		call quicksort(tout,isort,1,nt,near)
		rmed=tout((1+nt)/2)
		dev=0
		do n=1,nt
			dev=dev+abs(tout(n)-rmed)
		enddo
		dev=dev/nt
		devout(jp)=dev
		rmedout(jp)=rmed
		do n=1,nt
			outtr(ncp(n))=(outtr(ncp(n))-rmed)/dev
		enddo
	enddo !jp
	end
c
c	-------------------------------------------------------
	subroutine locateouts(clts,nearts,nclass,labelts,
     &	outts,jests,ncpts,ndpoints,jprox,npcase,kpr,kprox,
     &  near,jbt,touts,iout)
c
    	real outts(nearts),touts(nearts)
	integer nearts, nclass, labelts, near, jbt, n, jb, k, j, jp,
     &	jpclass, nt, les
c
	integer ncpts(nearts),clts(nearts),jests(nearts),
     &	ndpoints(nearts,jbt),iout(nearts),
     &  jprox(near,jbt),npcase(near,jbt),
     &	kpr(near),kprox(near,jbt),
     &	less(100)
c
	real rmed, amed, dev, adev, rqq
c
	if(labelts.eq.0) jpclass=0
	if(labelts.eq.1) then
		jpclass=nclass
		do n=1,nearts
			jests(n)=clts(n)
		enddo
	endif
	do n=1,nearts
		call zerv(kpr,near)
		do jb=1,jbt
			k=kprox(ndpoints(n,jb),jb)
			do j=jprox(k,jb),jprox(k+1,jb)-1
				kpr(npcase(j,jb))=kpr(npcase(j,jb))+1
			enddo
		enddo
		rqq=0
		do k=1,near
			if(kpr(k).gt.0) rqq=rqq+(kpr(k))**2
		enddo
		outts(n)=near/rqq
	enddo !n
	call zerv(less,100)
	les=0
	do jp=1,jpclass
		nt=0
		do n=1,nearts
			if(jests(n).eq.jp) then
				nt=nt+1
				touts(nt)=outts(n)
				ncpts(nt)=n
			endif
		enddo
		if(nt.gt.0.and.nt.le.4) then
			less(jp)=1
			les=1
			goto 48
		endif
		call quicksort(touts,iout,1,nt,nearts)
		rmed=touts((1+nt)/2)
		dev=0
		do n=1,nt
			dev=dev+abs(touts(n)-rmed)
		enddo
		dev=dev/nt
		do n=1,nt
			outts(ncpts(n))=(outts(ncpts(n))-rmed)/dev
		enddo
48		continue
	enddo !jp
	if(les.eq.1) then
		call quicksort(outts,iout,1,nearts,nearts)
		amed=outts(nint(.5*nearts))
		adev=0
		do n=1,nearts
			adev=adev+abs(outts(n)-amed)
		enddo
		adev=adev/nearts
		do jp=1,jpclass
			if (less(jp).eq.1) then
				do n=1,nt
					outts(ncpts(n))=(outts(ncpts(n))-amed)/adev
				enddo
			endif
		enddo
	endif
	end
c
c	-------------------------------------------------------
	subroutine myscale(loz,prox,xsc,y,u,near,nscale,red,nrnn,
     &	ee,ev,dl,kend)
c
	double precision prox(near,nrnn),y(near),u(near),dl(nscale),
     &	xsc(near,nscale),red(near),ee(near),ev(near,nscale),bl(10)
	integer loz(near,nrnn),kend(near)
c
	integer near, nscale, nrnn, j, i, it, n, jit, k
	double precision dotd, y2, sred, eu, ru, ra, ynorm
	do j=1,near
		ee(j)=dble(1.)
	enddo
c
	do j=1,near
		red(j)=0
		do i=1,nrnn
			red(j)=red(j)+prox(j,i)
		enddo
		red(j)=red(j)/near
	enddo
	sred=dotd(ee,red,near)
	sred=sred/near
	do it=1,nscale
		do n=1,near
			if(mod(n,2).eq.0) then
				y(n)=1
			else
				y(n)=-1
			endif
		enddo
		do jit=1,1000
			y2=dotd(y,y,near)
			y2=dsqrt(y2)
			do n=1,near
				u(n)=y(n)/y2
			enddo
			do n=1,near
				y(n)=0
				do k=1,kend(n)
					y(n)=y(n)+prox(n,k)*u(loz(n,k))
				enddo
			enddo
			eu=dotd(ee,u,near)
			ru=dotd(red,u,near)
			do n=1,near
				y(n)=y(n)-(red(n)-sred)*eu-ru
				y(n)=.5*y(n)
			enddo
			if(it.gt.1) then
				do j=1,it-1
					bl(j)=0
					do n=1,near
						bl(j)=bl(j)+ev(n,j)*u(n)
					enddo
					do n=1,near
						y(n)=y(n)-bl(j)*dl(j)*ev(n,j)
					enddo
				enddo
			endif
			ra=dotd(y,u,near)
			ynorm=0
			do n=1,near
				ynorm=ynorm+(y(n)-ra*u(n))**2
			enddo
			!write(*,*) jit,ynorm
			if(ynorm.lt.ra*1.0e-7)then
				do n=1,near
					xsc(n,it)=dsqrt(ra)*u(n)
					ev(n,it)=u(n)
				enddo
				dl(it)=ra
				goto 101
			endif
		enddo
101		continue
	enddo !nn
	end
c
c	-------------------------------------------------------
	subroutine fixfast(x,v,ncase,mdim,nsample,cat,code,
     &	nrcat,maxcat,fill,cl,nclass,kf,ff,jl,krcat)
c
c	input:
	real x(mdim,nsample),code
	integer cat(mdim),cl(nsample)
	integer mdim,nsample,maxcat,nclass
c	output:
	real fill(mdim,nclass)
c	work:
	real v(nsample),ff(nclass)
	integer nrcat(maxcat),krcat(maxcat)
	integer ncase(nsample),kf(mdim,nclass),jl(nclass)
c	local:
	integer m, jc, nt, jt, ii, jj, lcat, jjcat, n, j
	integer nmax, jmax, jjmax
c
	call zerm(kf,mdim,nclass)
	do m=1,mdim
	   if(cat(m).eq.1) then
		do jc=1,nclass
			nt=0
			do n=1,nsample
				if(abs(x(m,n)-code).gt.8.232D-11.and.cl(n).eq.jc) then
					nt=nt+1
					v(nt)=x(m,n)
				endif
			enddo
			call quicksort(v,ncase,1,nt,nsample)
			if(nt.gt.0) then
				fill(m,jc)=v((nt+1)/2)
			else
				kf(m,jc)=1
			endif
		enddo !jc
		jt=0
		do jj=1,nclass
			if(kf(m,jj).eq.0) then
				jt=jt+1
				ff(jt)=fill(m,jj)
			endif
		enddo!jj
		call quicksort(ff,jl,1,jt,nclass)
		do ii=1,nclass
			if(kf(m,ii).eq.1) fill(m,ii)=ff(jt/2)
		enddo!ii
	   endif !cat
	enddo !m
	do m=1,mdim
	   call zerv(krcat,maxcat)
	   if (cat(m).gt.1) then
		lcat=cat(m)
		call zerv(nrcat,maxcat)
		do jc=1,nclass
		do n=1,nsample
		if(abs(x(m,n)-code).gt.8.232D-11.and.cl(n).eq.jc) then
			j=nint(x(m,n))
			nrcat(j)=nrcat(j)+1
		endif
		enddo
		nmax=0
		jmax=1
		do j=1,lcat
		if(nrcat(j).gt.nmax) then
			nmax=nrcat(j)
			krcat(j)=krcat(j)+nrcat(j)
			jmax=j
		endif
		enddo
		if(nmax.ge.1) then
			fill(m,jc)=real(jmax)
		else
			kf(m,jc)=1
		endif
		enddo !jc
		jjcat=0    
		do j=1,lcat
			if(krcat(j).gt.jjcat) then
				jjcat=krcat(j)
				jjmax=j
			endif
		enddo
		do jc=1,nclass
			if(kf(m,jc).eq.1) fill(m,jc)=jjmax
		enddo
	   endif !cat
	enddo !m
	end
c
c	-------------------------------------------------------
	subroutine xfill(x,nsample,mdim,fill,cl,nclass,code)
c	
c	input:
	real code,fill(mdim,nclass)
	integer cl(nsample)
	integer mdim,nsample,nclass
c	output:
	real x(mdim,nsample)
c	local:
	integer n, m, jc
	do n=1,nsample
		do m=1,mdim
			do jc=1,nclass
				if(abs(x(m,n)-code).lt.8.232D-11.and.cl(n).eq.jc) 
     &				x(m,n)=fill(m,jc)
			enddo
		enddo !jc
	enddo !m
	end
c
c	-------------------------------------------------------
	subroutine roughfix(x,v,mdim,nsample,near,
     &	missing,cat,code,nrcat,maxcat,fill)
c
	real x(mdim,nsample), v(nsample), fill(mdim), code
	integer missing(mdim,nsample),cat(mdim),nrcat(maxcat)
	integer mdim, nsample, maxcat, near
	integer m, n, nt, j, jmax, lcat, nmax
	real rmed
	do n=1,near
		do m=1,mdim
			if(abs(code-x(m,n)).lt.8.232D-11) then
				missing(m,n)=1
			else 
				missing(m,n)=0
			endif
		enddo
	enddo
	do m=1,mdim
		if(cat(m).eq.1) then
c			continuous variable
			nt=0
			do n=1,nsample
				if(abs(x(m,n)-code).ge.8.232D-11) then
					nt=nt+1
					v(nt)=x(m,n)
				endif
			enddo
			if(nt.gt.0) then
				rmed=v((nt+1)/2)
				fill(m)=rmed
			else
				rmed=0
			endif
			do n=1,nsample
				if(abs(x(m,n)-code).lt.8.232D-11) x(m,n)=rmed
			enddo
		endif
		if (cat(m).gt.1) then
c			categorical variable
			lcat=cat(m)
			call zerv(nrcat,maxcat)
			do n=1,nsample
				if(abs(x(m,n)-code).ge.8.232D-11) then
					j=nint(x(m,n))
					nrcat(j)=nrcat(j)+1
				endif
			enddo
			nmax=0
			jmax=1
			do j=1,lcat
				if(nrcat(j).gt.nmax) then
					nmax=nrcat(j)
					jmax=j
				endif
			enddo
			fill(m)=real(jmax)
			do n=1,nsample
				if(abs(x(m,n)-code).lt.8.232D-11) x(m,n)=real(jmax)
			enddo
		endif
	enddo !m
	end
c
c	-------------------------------------------------------
	subroutine fixtestunlab(xxts,xts,mdim,nclass,ntest,fill,code)
c
	real xts(mdim,ntest),xxts(mdim,nclass*ntest),fill(mdim,nclass),code
	integer mdim, nclass, ntest
	integer n, m, j, nn
c	unlabelled test set - get a fill value for each possible class
	do n=1,ntest
	do j=1,nclass
		nn=(n-1)*nclass+j
		do m=1,mdim
		if(abs(xts(m,n)-code).lt.8.232D-11) then
			xxts(m,nn)=fill(m,j)
		else
			xxts(m,nn)=xts(m,n)
		endif
		enddo
	enddo
	enddo
	end
c
c	-------------------------------------------------------
	subroutine fixtestlab(xts,mdim,nclass,ntest,fill,code,clts)
c
	real xts(mdim,ntest),fill(mdim,nclass),code
	integer mdim, nclass, ntest, clts(ntest)
	integer n, m
c	labelled test set - use the training set fill value (for this var/class)
	do n=1,ntest
		do m=1,mdim
			if(abs(xts(m,n)-code).lt.8.232D-11) then
				xts(m,n)=fill(m,clts(n))
			endif
		enddo
	enddo
	end
c
c	-------------------------------------------------------
	subroutine reducemult(wr,nclass,ntest,qxxts,jests,qts,
     &	xxts,xts,mdim,imax)
c
	real wr(nclass),qxxts(nclass,nclass*ntest),
     &	qts(nclass,ntest),xxts(mdim,nclass*ntest),xts(mdim,ntest)
   	integer nclass,ntest,mdim,jests(ntest),imax(ntest)
	integer n, j, i, m, jj, nn
	real qmax 
c
	do n=1,ntest
		call zervr(wr,nclass)
		do j=1,nclass
			nn=(n-1)*(nclass)+j
			do i=1,nclass
				wr(i)=qxxts(i,nn)
			enddo !i
		enddo !i
		qmax=0
		do i=1,nclass
			if(wr(i).ge.qmax) then
				imax(n)=i
				qmax=wr(i)
			endif
		enddo
		jests(n)=imax(n)
		do m=1,mdim
			xts(m,n)=xxts(m,(n-1)*nclass+imax(n))
		enddo
		do jj=1,nclass
			qts(jj,n)=wr(jj)
		enddo !jj
	enddo !n
	end
c
c	-------------------------------------------------------
	subroutine impute(x,prox,near,mdim,nclass,
     &	maxcat,votecat,cat,nrnn,loz,cl,kend,
     &	missing)
c
	real x(mdim,near),votecat(maxcat)
	double precision  prox(near,nrnn)
	integer near, mdim, nclass, maxcat, nrnn
	integer cat(mdim),loz(near,nrnn),cl(near),
     &	kend(near),missing(mdim,near)
	integer i, j, jmax, m, jc, n, k
	real sx, dt, rmax
c
	do m=1,mdim
		if(cat(m).eq.1) then
			do jc=1,nclass
				do n=1,near
					if(missing(m,n).eq.1.and.cl(n).eq.jc) then
						sx=0
						dt=0
						do k=1,min0(nrnn,kend(n))
							if (missing(m,loz(n,k)).ne.1.and.
     &							cl(loz(n,k)).eq.jc) then
								sx=sx+real(prox(n,k))*x(m,loz(n,k))
								dt=dt+real(prox(n,k))
							endif
						enddo
						if(dt.gt.0) x(m,n)=sx/dt
					endif
				enddo !n
			enddo !jc
		endif
	enddo !m
	do m=1,mdim
		if(cat(m).gt.1) then
			do jc=1,nclass
			do n=1,near
				if(missing(m,n).eq.1.and.cl(n).eq.jc) then
					call zervr(votecat,maxcat)
					do k=1,min0(nrnn,kend(n))
						if (missing(m,loz(n,k)).ne.1.and.
     &							cl(loz(n,k)).eq.jc) then
							j=nint(x(m,loz(n,k)))
							votecat(j)=votecat(j)+real(prox(n,k))
						endif
					enddo !k
					rmax=-1
					do i=1,cat(m)
						if(votecat(i).gt.rmax) then
							rmax=votecat(i)
							jmax=i
						endif
					enddo
					x(m,n)=real(jmax)
				endif
			enddo !n
			enddo !jc
		endif
	enddo !m
	end
c
c	-------------------------------------------------------
	subroutine quicksort(v,iperm,ii,jj,kk)
c
c	puts into iperm the permutation vector which sorts v into
c	increasing order.  only elementest from ii to jj are considered.
c	array iu(k) and array il(k) permit sorting up to 2**(k+1)-1 elements
c
c	this is a modification of acm algorithm #347 by r. c. singleton,
c	which is a modified hoare quicksort.
c
	real v(kk), vt, vtt
	integer t,tt,iperm(kk),iu(32),il(32)
	integer ii, jj, kk, m, i, j, k, ij, l
c
	m=1
	i=ii
	j=jj
 10   	if (i.ge.j) go to 80
 20   	k=i
	ij=(j+i)/2
	t=iperm(ij)
	vt=v(ij)
	if (v(i).le.vt) go to 30
	iperm(ij)=iperm(i)
	iperm(i)=t
	t=iperm(ij)
	v(ij)=v(i)
	v(i)=vt
	vt=v(ij)
 30   	l=j
	if (v(j).ge.vt) go to 50
	iperm(ij)=iperm(j)
	iperm(j)=t
	t=iperm(ij)
	v(ij)=v(j)
	v(j)=vt
	vt=v(ij)
	if (v(i).le.vt) go to 50
	iperm(ij)=iperm(i)
	iperm(i)=t
	t=iperm(ij)
	v(ij)=v(i)
	v(i)=vt
	vt=v(ij)
	go to 50
 40   	iperm(l)=iperm(k)
	iperm(k)=tt
	v(l)=v(k)
	v(k)=vtt
 50   	l=l-1
	if (v(l).gt.vt) go to 50
	tt=iperm(l)
	vtt=v(l)
 60   	k=k+1
	if (v(k).lt.vt) go to 60
	if (k.le.l) go to 40
	if (l-i.le.j-k) go to 70
	il(m)=i
	iu(m)=l
	i=k
	m=m+1
	go to 90
 70   	il(m)=k
	iu(m)=j
	j=l
	m=m+1
	go to 90
 80   	m=m-1
	if (m.eq.0) return
	i=il(m)
	j=iu(m)
 90   	if (j-i.gt.10) go to 20
	if (i.eq.ii) go to 10
	i=i-1
 100  	i=i+1
	if (i.eq.j) go to 80
	t=iperm(i+1)
	vt=v(i+1)
	if (v(i).le.vt) go to 100
	k=i
 110  	iperm(k+1)=iperm(k)
	v(k+1)=v(k)
	k=k-1
	if (vt.lt.v(k)) go to 110
	iperm(k+1)=t
	v(k+1)=vt
	go to 100
	end
c
c	-------------------------------------------------------
	subroutine quicksortint(v,iperm,ii,jj,kk)
c
c	puts into iperm the permutation vector which sorts v into
c	increasing order.  only elementest from ii to jj are considered.
c	array iu(k) and array il(k) permit sorting up to 2**(k+1)-1 elements
c
c	v is also sorted
c
c	this is a modification of acm algorithm #347 by r. c. singleton,
c	which is a modified hoare quicksort.
c
	integer v(kk),vt,vtt
	integer t,tt,iperm(kk),iu(32),il(32)
	integer ii, jj, m, i, j, k, ij, kk, l
	m=1
	i=ii
	j=jj
 10   	if (i.ge.j) go to 80
 20   	k=i
	ij=(j+i)/2
	if(ij.gt.kk.or.ij.lt.1) write(*,*) ij,ii,kk
	t=iperm(ij)
	vt=v(ij)
	if (v(i).le.vt) go to 30
	iperm(ij)=iperm(i)
	iperm(i)=t
	t=iperm(ij)
	v(ij)=v(i)
	v(i)=vt
	vt=v(ij)
 30   	l=j
	if (v(j).ge.vt) go to 50
	iperm(ij)=iperm(j)
	iperm(j)=t
	t=iperm(ij)
	v(ij)=v(j)
	v(j)=vt
	vt=v(ij)
	if (v(i).le.vt) go to 50
	iperm(ij)=iperm(i)
	iperm(i)=t
	t=iperm(ij)
	v(ij)=v(i)
	v(i)=vt
	vt=v(ij)
	go to 50
 40   	iperm(l)=iperm(k)
	iperm(k)=tt
	v(l)=v(k)
	v(k)=vtt
 50   	l=l-1
	if (v(l).gt.vt) go to 50
	tt=iperm(l)
	vtt=v(l)
 60   	k=k+1
	if (v(k).lt.vt) go to 60
	if (k.le.l) go to 40
	if (l-i.le.j-k) go to 70
	il(m)=i
	iu(m)=l
	i=k
	m=m+1
	go to 90
 70   	il(m)=k
	iu(m)=j
	j=l
	m=m+1
	go to 90
 80   	m=m-1
	if (m.eq.0) return
	i=il(m)
	j=iu(m)
90	if (j-i.gt.10) go to 20
	if (i.eq.ii) go to 10
	i=i-1
100	i=i+1
	if (i.eq.j) go to 80
	t=iperm(i+1)
	vt=v(i+1)
	if (v(i).le.vt) go to 100
	k=i
110	iperm(k+1)=iperm(k)
	v(k+1)=v(k)
	k=k-1
	if (vt.lt.v(k)) go to 110
	iperm(k+1)=t
	v(k+1)=vt
	go to 100
	end
c
c	-------------------------------------------------------
	subroutine unpack(l,npack,icat)
c
	integer icat(32),npack,l,j,n,k
	if(l.gt.32) then
		write(*,*) "error in unpack, l=",l
		stop
	endif
	do j=1,32
		icat(j)=0
	enddo
	n=npack
	icat(1)=mod(n,2)
	do k=2,l
		n=(n-icat(k-1))/2
		icat(k)=mod(n,2)
	enddo
	end
c
c	-------------------------------------------------------
	subroutine zerv(ix,m1)
c
	integer ix(m1), m1, n
	do n=1,m1
		ix(n)=0
	enddo
	end
c
c	-------------------------------------------------------
	subroutine zervr(rx,m1)
c
	real rx(m1)
	integer m1,n
	do n=1,m1
		rx(n)=0
	enddo
	end
c
c	-------------------------------------------------------
	subroutine zerm(mx,m1,m2)
c
	integer mx(m1,m2),m1,m2,i,j
	do j=1,m2
		do i=1,m1
			mx(i,j)=0
		enddo
	enddo
	end
c
c	-------------------------------------------------------
	subroutine zermr(rx,m1,m2)
c
	real rx(m1,m2)
	integer m1,m2,i,j
	do j=1,m2
		do i=1,m1
			rx(i,j)=0
		enddo
	enddo
	end
c
c	-------------------------------------------------------
	subroutine zermd(rx,m1,m2)
c
	double precision rx(m1,m2)
	integer m1,m2,i,j
	do j=1,m2
		do i=1,m1
			rx(i,j)=0
		enddo
	enddo
	end
c
c	-------------------------------------------------------
	subroutine eqm(j,k,m,n,mt,msm)
c
c	copies rows msm(1),...,msm(mdimt) of k into the same rows of j
	integer j(m,n),k(m,n),m,n,mt,msm(mt)
	integer n1,k1,m1
	do n1=1,n
		do k1=1,mt
			m1=msm(k1)
			j(m1,n1)=k(m1,n1)
		enddo
	enddo
	end
c
c	-------------------------------------------------------
	real function erfcc(x)
c
	real x,t,z
	z=abs(x)/1.41421356
	t=1./(1.+0.5*z)
	erfcc=t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+t*
     *	(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+t*
     *	(1.48851587+t*(-.82215223+t*.17087277)))))))))
	erfcc=erfcc/2
	if (x.lt.0.) erfcc=2.-erfcc
	return
	end
c
c	-------------------------------------------------------
	real function rrand(i)
c
	double precision dseed,u
	integer i
	save dseed
	data dseed /17395/
	call lrnd(dseed,u)
	rrand=real(u)
	end
c
c	-------------------------------------------------------
	subroutine lrnd(dseed,u)
c
	double precision dseed, d31m1,u
	data d31m1 /2147483647/
	dseed=dmod(16087*dseed,d31m1)
	u=dseed/d31m1
	return
	end
c
c	-------------------------------------------------------
	real function rnorm(i)
c
	integer i
	real rrand
	rnorm=sqrt(-2*log(rrand(1)))*cos(6.283185*rrand(1))
	end 
c
c	-------------------------------------------------------
	double precision function dotd(u,v,ns)
c
c	computes the inner product
c	input:
	double precision u(ns),v(ns)
	integer ns
c	local:
	integer n
	dotd=dble(0.0)
	do n=1,ns
		dotd=dotd+u(n)*v(n)
	enddo
	end
