program randomforestsV1

C     Copyright (C)  2001 Adele Cutler and Leo Breiman
C     This program is free software; you can redistribute it and/or modify
C     it under the terms of the GNU General Public License as published by
C     the Free Software Foundation; either version 2 of the License, or (at
C     your option) any later version.
C     This program is distributed in the hope that it will be useful, but
C     WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
C     or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
C     for more details.                   
        

C       SEE THE GUIDE TO USING RANDOM FORESTS FOR INSTRUCTIONS 

 
c       set all parameters except those in last line.  

c       line #1 information about the data
c       line #2 setting up the run
c       line #3 options
c       line #4 priors and costs
c       line #5 ooutput controlS
c       line #6 memory parameters (automatically set)

        parameter(
     1  mdim=9,nsample=200,nclass=6,maxcat=1,ntest=14,ltest=0,
     2  mtry=3,jbt=100,look=1,modevar=0,
     3  imp=0,igini=0,iprox=0,iden=0,noutlier=0,
     4  ipi=1,icost=1,
     5  infow=0,impw=0,iginiw=0,iproxw=0,idenw=0,noutlierw=0,ntestw=1,
     6	nrnodes=2*(nsample)+1,mimp=imp*(mdim-1)+1,nimp=imp*(nsample-1)+1,
     6	near=iprox*(nsample-1)+1, nlier=noutlier*(nsample-1)+1)
     
        
c       MAIN PROGRAM

c       DIMENSIONING OF ARRAYS


        real    x(mdim,nsample),pi(nclass),misscost(nclass,nclass),
     1       sm(mdim),rcat(maxcat,nrnodes),v(nsample),ut(nsample),
     1       xt(nsample),xb(mdim,nsample),tp(nsample),tx(nsample),
     1       errimp(mimp),pid(nclass),pidsq(nclass),upper(nrnodes),
     1       rmargin(nsample),tgini(mdim),outlier(nlier),
     1       q(nclass,nsample),prox(near,near),pr(mdim),tsuts(ntest),
     1       oulierclass(nclass,nlier),uts(nsample),xts(mdim,ntest)
     1      
                
                
        
        integer cl(nsample),jdex(nsample),jerr(nsample),mbest(nrnodes),
     1       treemap(2,nrnodes), classpop(nclass,nrnodes),jl(nsample),
     1       nodestatus(nrnodes),nodepop(nrnodes),jts(nsample),
     1       clb(nsample),nodeclass(nrnodes),nperm(nsample),
     1       parent(nrnodes),nc(nclass),cat(mdim),nout(nsample),
     1       jest(nsample),itc(maxcat),lcl(nsample),jin(nsample),
     1       jtr(nsample),ncounttr(nclass,nsample),lc(maxcat),
     1       missimp(mimp),isort(nsample),ncountimp(nclass,nimp,mimp),
     1       tclasspop(nclass),nodestart(nrnodes),tclassr(nclass),
     1       tclassl(nclass),ncase(nsample),tmpclass(nclass),
     1       tlsclasspop(nclass),iv(mdim),incl(mdim),
     1       nodex(nsample),clts(ntest),jtts(ntest),jmaxts(ntest),
     1       ncount(nclass,ntest),nodexts(ntest),mtab(nclass,nclass)
                
                        
                
        

        
c       initialize seed in random number generator

        do n=1,711
        zz=randn(1)
        end do

c       BEGIN MAIN
        
C       READ IN THE DATA+++++++++++++++++++++++++++++++++++++++
	nset=2
c       cl(k) is the class label, declared integer. x(j,k)  are the input variables
        if(nset.eq.1) then
        open(9,file='xxxx.train',status='old')
        do k=1,nsample
        read(9,*) cl(k),(x(j,k),j=1,mdim)
        end do
        close(9)
        
        if(ltest.gt.0) then
        open(7,file='xxxx.test',status='old')
        do k=1,ntest
        read(7,*) clts(k),(xts(j,k),j=1,mdim)
        end do
        close(7)
        end if
	end if
        
	if (nset.eq.2) then
	open(9,file='glass6',status='old')
	do k=1,200
	read(9,*) cl(k),extra,(x(j,k),j=1,mdim)
	end do
	do k=1,ntest
	read(9,*) clts(k),extra,(xts(j,k),j=1,mdim)
	end do
	close(9)
	end if
	
	
        
c       SET CATEGORICAL VALUES+++++++++++++++++++++++++++++++++++
        
        if(maxcat.eq.1) then
        do m=1,mdim
        cat(m)=1
        end do
        else
        do m=1,mdim
        !cat(m)=999
        end do
        end if
	do m=1,mdim
	!write(*,*) m,cat(m)
	end do

        call setprob(pr,sm,mdim,cat)
        
c       SET WHICH VARIABLES TO INCLUDE+++++++++++++++++++++++++++++++

        if(modvar.eq.0) then
        do m=1,mdim
        incl(m)=1
        end do
        end if

        if (modvar.eq.1) then
        do m=1,mdim
        incl(m)=0
        end do
        do m=1,13
        incl(m)=1
        end do
        end if
        
c       SET UP DATA FOR DENSITY OR OUTLIERS++++++++++++++++++++++++++++++

        if(iden.eq.1.or.noutlier.eq.1) then
        nhalf=nsample/2
        do n=1,nhalf
        cl(n)=1
        end do
        do n=nhalf+1,nsample
        cl(n)=2
        end do
        do n=nhalf+1,nsample
        do m=1,mdim
        k=int(randn(1)*nhalf)+1
        x(m,n)=x(m,k)
        end do
        end do
        end if
        
c       initialize+++++++++++++++++++++++++++++++++++++++++++++++++++++++

        averrc=0
        
        call zerm(ncounttr,nclass,nsample)
        call zerv(nout,nsample)
        call zermr(prox,near,near)
        call zervr(tgini,mdim)
        
        if(imp.eq.1) then
        do j=1,nclass
        do n=1,nsample
        do m=1,mdim
        ncountimp(j,n,m)=0
        end do
        end do
        end do
        end if
        
        
        
c       BEGIN CONSTRUCTING TREES=========================================

        do jb=1,jbt
        
        
c       fill in omitted variables by noise

        if (modvar.eq.1) then
        do m=1,mdim
        if(incl(m).eq.0) then
        do n=1,nsample
        xt(n)=x(m,n)
        end do
        call perm1(nsample,nsample,xt)
        do n=1,nsample
        x(m,n)=xt(n)
        end do
        end if
        end do
        end if
        
        
        
c       set up bootstrap training set

51      continue
        call zerv(jin,nsample)
        call zerv(nc,nclass)
        do n=1,nsample
        k=int(randn(1)*nsample)+1
        jin(k)=1
        clb(n)=cl(k)
        nc(cl(k))=nc(cl(k))+1
        do m=1,mdim
        xb(m,n)=x(m,k)
        end do
        end do
        
        
        
        ncmax=0
        do j=1,nclass
        ncmax=max0(ncmax,nc(j))
        end do
        if(ncmax.eq.nsample) goto 51

c       build tree

        call prep(nsample,clb,nclass,ipi,pi,pid,pidsq,tclasspop,icost,
     1  misscost)
        
        call buildtree(xb,clb,mdim,nsample,nclass,pid,pidsq,treemap,
     1  jdex,upper,nodestatus,nodepop,nodestart,classpop,tclasspop,
     1  tclassl,tclassr,nrnodes,ncase,parent,tmpclass,ut,v,tgini,
     1  ndbigtree,xt,jl,mtry,mbest,cat,maxcat,itc,rcat,sm,nodeclass,lc,
     1  misscost,iv)
        
        
c       ++++++++++++++++++++++++++++++++        
        
c       compute oob stuff-test set error, estimated class prob, etc.

        call testreebag(x,nsample,mdim,treemap,nodestatus,nodeclass,
     1  nrnodes,ndbigtree,nclass,jtr,uts,upper,mbest,rcat,cat,maxcat,
     1  nodex)
        
        call oob(nsample,nclass,jin,cl,jtr,jerr,ncounttr,nout,errtr,
     1  errc,jout,rmargin,q,jest)
     
        
        averrc=((jb-1)*averrc+errc)/jb
        
c       do variable importance.++++++++++ 

        if(imp.eq.1) then
        do mr=1,mdim
        if(incl(mr).eq.1) then
        call permobmr(mr,x,tp,tx,jin,nsample,mdim)
        call testreebag(x,nsample,mdim,treemap,nodestatus,nodeclass,
     1  nrnodes,ndbigtree,nclass,jts,uts,upper,mbest,rcat,cat,maxcat,
     1  nodex)
        call contimp(mr,x,jin,jts,tx,ncountimp,nclass,nsample,mdim)
        end if
        end do !mr
        call finishimp(missimp,ncountimp,nout,cl,nclass,mdim,nsample,
     1  errimp,jout)
        end if !imp
        
        
c       evaluate proximities+++++++++++++

        if(iprox.eq.1) then
        call ndnear(nsample,jin,nodex,prox)
        end if
        
c       work with test set++++++++++++

        if(ltest.gt.0) then
        
        call testreebag(xts,ntest,mdim,treemap,nodestatus,nodeclass,
     1  nrnodes,ndbigtree,nclass,jtts,tsuts,upper,mbest,rcat,cat,maxcat,
     1  nodexts)
        
	call comptserr(ncount,jtts,jmaxts,clts,ntest,nclass,errts,ltest)
        end if

        

c       write output on run+++++++++++++++

        if(mod(jb,look).eq.0) then
        if(ltest.le.1)then
        write(*,*) jb, 100*errtr
        end if
	if(ltest.eq.2) then
        write(*,*) jb,100*errts,100*errtr
        end if
        end if
        
        end do !jbt
        
c       END TREE CONSTRUCTION==========================================

c       OUTPUT  OUTPUT OUTPUT==========================================

        
        write(*,*) ""
        if(ltest.lt.2)then
        write(*,*) "final error rate %    ",100*errtr
        end if
	if (ltest.eq.2) then
        write(*,*) "final error rate %    ",100*errtr
        write(*,*) "final error test %    ",100*errts
	end if
        
        write(*,*) ""
        call zerm(mtab,nclass,nclass)
        do n=1,nsample
        mtab(cl(n),jest(n))=mtab(cl(n),jest(n))+1
        end do
        write(*,*) "          true class "
        write(*,*) ""
        write(*,'(20i6)') "", (i, i=1, nclass)
        write(*,*) ""
        do j=1,nclass
        write(*,'(20i6)')  j,(mtab(i,j), i=1,nclass)
        end do
        write(*,*) ""
        
        
C       SEND INFOUT DATA TO FILE++++++++++++++++++++++++++++++++++++++
        
        if(infow.eq.1) then
        open(7,file='info-xxxx', status='new')
        do n=1,nsample
        write(7,'(4i5,20f10.3)') n,jerr(n),cl(n),jest(n),rmargin(n),
     1  (q(j,n), j=1,nclass)
        end do
        cLose(7)
        end if
        
C       SEND PROXIMITY DATA TO FILE++++++++++++++++++++++++++++++++
        
        if(iprox.eq.1.or.noutlierw.eq.1) then
        do n=1,nsample
        do k=1,n
        prox(n,k)=(prox(n,k)+prox(k,n))/jbt
        if(k.eq.n) prox(n,k)=0
        end do
        end do
        do n=1,nsample
        do k=1,n
        prox(k,n)=prox(n,k)
        end do
        end do
        end if
        if(iproxw.eq.1) then 
        open(9, file='prox-xxxx',status='new')
        do n=1,nsample
        do k=1,nsample
        write(9,'(2i5,f10.3)') n,k,prox(n,k)
        end do
        end do
        close(9)
        end if
        
C       SEND DENSITY DATA TO FILE++++++++++++++++++++++++++++++++
        
        if(idenw.eq.1) then
        open(9,file='density-xxxx',status='new')
        do n=1,nhalf
        write(9,'(i5,3f10.3)') n, q(1,n),q(2,n),q(1,n)/(q(2,n)+1.0e-6)
        end do
        end if
        
C       SEND OUTLIER DATA TO FILE+++++++++++++++++++++++++++++
        
        
        if(noutlierw.eq.1) then
        open(9,file='outlier-xxxx',status='new')
        do n=1,nsample/2
        outlier(n)=0
        do k=1,nsample/2
        outlier(n)=outlier(n)+(prox(n,k))**2
        end do
        end do
        do n=1,nsample/2
        outlier(n)=1.0/outlier(n)
        end do
        do n=1,nsample/2
        write(9,'(i5,f10.2)') n,outlier(n)
        end do
        close(9)
        end if
        
C       SEND IMP DATA TO FILE++++++++++++++++++++++++++++++++++
        
        if (impw.eq.1) then
        open(9,file='importance-xxxx',status='new')
        do m=1,mdim
        if(incl(m).eq.1) then
        write(9,'(i5,f10.2)')m,100*(errimp(m)-errtr)/errtr
        end if
        end do
        close(9)
        end if
        
C       SEND GINI DATA TO FILE++++++++++++++++++++++++++++++
        
        if(iginiw.eq.1) then
	avg=0
	sqg=0
	do m=1,mdim
	sqg=sqg+(tgini(m)-avg)**2
	avg=((m-1)*avg+tgini(m))/m
	end do
	sd=sqrt(sqg/mdim)
	do m=1,mdim
	tgini(m)=(tgini(m)-avg)/sd
	if(tgini(m).le.0) tgini(m)=0
	end do
        open(9,file='gini-xxxx',status='new')
        do m=1,mdim
        write(9,'(i5,f10.2)') m,100*tgini(m)
        end do
        close(9)
        end if
        
c       SEND TEST SET DATA TO FILE++++++++++++++++++++++++++++

        If (ntestw.gt.0) then
        !open(9,file='testout-xxxx',status='new')
        do n=1,ntest
        write(*,'(3i5,20f10.2)') n,clts(n),jmaxts(n),
     1  (real(ncount(j,n))/jbt, j=1,nclass)
        end do
        end if
        
        
        
        end !MAIN
                
        
c       END MAIN

        subroutine setprob(pr,sm,mdim,cat)
        real pr(mdim),sm(mdim)
        integer cat(mdim)
        mc=0
        do m=1,mdim
        if (cat(m).le.2) then
        mc=mc+1
        pr(m)=1
        else
        mc=mc+cat(m)-1
        pr(m)=cat(m)-1
        end if
        end do
        
        do m=1,mdim
        pr(m)=pr(m)/mc
        end do
        sm(1)=pr(1)
        do m=2,mdim
        sm(m)=sm(m-1)+pr(m)
        end do
        end
        

c       SUBROUTINE BUILDTREE
        
        subroutine buildtree(x,cl,mdim,nsample,nclass,pid,pidsq,treemap,
     1       jdex,upper,nodestatus,nodepop,nodestart,classpop,tclasspop,
     1       tclassl,tclassr,nrnodes,ncase,parent,tmpclass,ut,v,tgini,
     1       ndbigtree,xt,jl,mtry,mbest,cat,maxcat,itc,rcat,sm,
     1       nodeclass,lc,misscost,iv)
        
        integer cl(nsample),treemap(2,nrnodes),tclasspop(nclass),
     1       parent(nrnodes),nodestatus(nrnodes),tmpclass(nclass),
     1       jl(nsample),nodepop(nrnodes),nodestart(nrnodes),
     1       classpop(nclass,nrnodes),tclassl(nclass),tclassr(nclass),
     1       jdex(nsample),ncase(nsample),cat(mdim),itc(maxcat),
     1       lc(maxcat),mbest(nrnodes),nodeclass(nrnodes),iv(mdim)
        
        real  pid(nclass),pidsq(nclass),x(mdim,nsample),
     1       xt(nsample),upper(nrnodes),v(nsample),sm(mdim),tgini(mdim),
     1       ut(nsample),rcat(maxcat,nrnodes),misscost(nclass,nclass)
        
        
        
        call zerv(nodestatus,nrnodes)
        call zerv(nodestart,nrnodes)
        call zerv(nodepop, nrnodes)
        call zerm(classpop,nclass,nrnodes)
        
        re=.5
        iprint=0
        
        do 20 j=1,nclass
        classpop(j,1)=tclasspop(j)
20      continue

        do n=1,nsample
        ut(n)=0
        jdex(n)=n
        end do

        ncur=1
        nodestart(1)=1
        nodepop(1)=nsample
        nodestatus(1)=2
        
c       start main loop

        do 30 kbuild=1,nrnodes
c       write(*,*) kbuild,(classpop(j,kbuild), j=1,nclass)
        if (kbuild.gt.ncur) goto 50
        if (nodestatus(kbuild).ne.2) goto 30
        
c               initialize for next call to findbestsplit

        ndstart=nodestart(kbuild)
        ndend=ndstart+nodepop(kbuild)-1
        do j=1,nclass
        tclasspop(j)=classpop(j,kbuild)
        end do
        
        iprint=0
        if(iprint.eq.1) then
        write(*,*) kbuild
        write(*,*) (classpop(j,kbuild), j=1,nclass)
        write(*,*) ndstart,ndend
        write(*,*) ""
        end if
        

        

        
        call findbestsplit(x,xt,ut,jdex,jl,cl,mdim,nsample,nclass,cat,
     1  maxcat,pid,pidsq,ndstart,ndend,tclasspop,tclassl,tclassr,msplit,
     1  decsplit,ubest,ncase,tmpclass,ndendl,jstat,v,mtry,itc,sm,kbuild,
     1  lc,nbestvar,iv)
        
        
        if (jstat.eq.1) then
        nodestatus(kbuild)=-1
        go to 30
        else
        mv=nbestvar
        tgini(mv)=decsplit+tgini(mv)
        mbest(kbuild)=mv
        if(cat(mv).gt.1) then
                do i=1,cat(mv)
                rcat(i,kbuild)=real(itc(i))
                end do
                end if
        upper(kbuild)=ubest
        end if
        
        
        
c       leftnode no.= ncur+1, rightnode no. = ncur+2.

        nodepop(ncur+1)=ndendl-ndstart+1
        nodepop(ncur+2)=ndend-ndendl
        nodestart(ncur+1)=ndstart
        nodestart(ncur+2)=ndendl+1

c               find class populations in both nodes
        
        do  n=ndstart,ndendl
        j=cl(jdex(n))
        classpop(j,ncur+1)=classpop(j,ncur+1)+1
        end do

        do n=ndendl+1,ndend
        j=cl(jdex(n))
        classpop(j,ncur+2)=classpop(j,ncur+2)+1
        end do
        
        
        

c       check on nodestatus
        
        nodestatus(ncur+1)=2
        nodestatus(ncur+2)=2
        if (nodepop(ncur+1).le.1) nodestatus(ncur+1)=-1
        if (nodepop(ncur+2).le.1) nodestatus(ncur+2)=-1
        do 80 j=1,nclass
        if (classpop(j,ncur+1).eq.nodepop(ncur+1)) nodestatus(ncur+1)=-1
        if (classpop(j,ncur+2).eq.nodepop(ncur+2)) nodestatus(ncur+2)=-1
80      continue
        
        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
        end do
        
        
        
        do ncost=1,ndbigtree
        pmin=1e20
        do iclass=1,nclass
        p=0
        do jclass=1,nclass
        p=p+misscost(iclass,jclass)*pid(jclass)*classpop(jclass,ncost)
        end do
        if (p.lt.pmin) then
        nodeclass(ncost)=iclass
        pmin=p
        endif
        end do
        end do
        

        end


c       SUBROUTINE FINDBESTSPLIT

        subroutine findbestsplit(x,xt,ut,jdex,jl,cl,mdim,nsample,nclass,
     1  cat,maxcat,pid,pidsq,ndstart,ndend,tclasspop,tclassl,tclassr,
     1  msplit,decsplit,ubest,ncase,tmpclass,ndendl,jstat,v,mtry,itc,sm,
     1  kbuild,lc,nbestvar,iv)
     1  
     
        integer cl(nsample),tclasspop(nclass),tmpclass(nclass),iv(mdim),
     1       tclassl(nclass),tclassr(nclass),ncase(nsample),cat(mdim),
     1       jdex(nsample),jl(nsample),lc(maxcat),itc(maxcat)
        
                
        real pid(nclass),pidsq(nclass),x(mdim,nsample),ut(nsample),
     1       xt(nsample),v(nsample),sm(mdim)        
        
c               compute initial values of numerator and denominator of Gini
        pno=0
        pdo=0
        do 10 j=1,nclass
        pno=pno+pidsq(j)*tclasspop(j)*tclasspop(j)
        pdo=pdo+pid(j)*tclasspop(j)
10      continue
        crit0=pno/pdo
        jstat=0
        critmax=-1.0e10
        critvar=-1.0e10
        non=0
	
	
                
                
                
100             call zerv(iv,mdim)              
                do it=1,mtry

200             if(maxcat.eq.1) then
                kv=int(mdim*randn(1))+1
                if(kv.gt.mdim) kv=mdim
                else            
                kv=nselect(mdim,sm)
		end if
                if(iv(kv).ne.0) goto 200
                iv(kv)=1
                if(cat(kv).gt.1) then
                jj=cat(kv)
                call catran(lc,jj,maxcat)
                end if
                
                
                
                do n=ndstart,ndend
                xt(n)=0
                if(cat(kv).eq.1) then
                xt(n)=x(kv,jdex(n))
                else
                xt(n)=lc(nint(x(kv,jdex(n))))
		end if
                jl(n)=cl(jdex(n))
                end do
                
                
                
                do n=ndstart,ndend
                v(n)=xt(n)
                end do
                do n=1,nsample
                ncase(n)=n
                end do
                call quicksort (v,ncase,ndstart,ndend,nsample)
                if(v(ndstart).ge.v(ndend)) then
                non=non+1
                if(non.gt.1000) then
                jstat=1
                return
                end if
                goto 100
                end if
                
c               ncase(n)=case number of v nth from bottom


                rrn=pno
                rrd=pdo
                rln=0
                rld=0
                call zerv(tclassl,nclass)
                do j=1,nclass
                tclassr(j)=tclasspop(j)
                end do

                do nsp=ndstart,ndend-1
                kc=jl(ncase(nsp))
                rln=rln+pidsq(kc)*(1+2*tclassl(kc))
                rrn=rrn+pidsq(kc)*(1-2*tclassr(kc))
                rld=rld+pid(kc)
                rrd=rrd-pid(kc)
                crit=(rln/rld)+(rrn/rrd)
                
                
                if (crit.gt.critvar) then
                        if (v(nsp).lt.v(nsp+1)) then
                                ubestt=v(nsp)+.5*(v(nsp+1)-v(nsp))
                                critvar=crit
                                nbestt=nsp
                                endif
                end if
                
                
                tclassl(kc)=tclassl(kc)+1
                tclassr(kc)=tclassr(kc)-1
                end do !nsp
                if(critvar.gt.critmax) then
                critmax=critvar
                ubest=ubestt
                nbest=nbestt
                nbestvar=kv
                if(cat(kv).gt.1) then
                do i=1,cat(kv)
                itc(i)=lc(i)
                end do
                end if
                do n=ndstart,ndend
                ut(n)=xt(n)
                end do
                end if
                
                end do !mtry
                
                
                nl=ndstart-1
                do nsp=ndstart,ndend
                if(ut(nsp).le.ubest) then
                nl=nl+1
                ncase(nl)=jdex(nsp)
                end if
                end do
                ndendl=max(nl,ndstart)
                nr=ndendl
                do nsp=ndstart,ndend
                if(ut(nsp).gt.ubest) then
                nr=nr+1
                if(nr.gt.nsample) goto 765
                ncase(nr)=jdex(nsp)
                end if
                end do
765             continue
                if(ndendl.eq.ndend) ndendl=ndend-1
                decsplit=critmax-crit0
                do n=ndstart,ndend
                jdex(n)=ncase(n)
                end do
                
                end
                

        


        
c       SMALL ASSOCIATED SUBROUTINES

        

        
        subroutine prep(nsample,cl,nclass,ipi,pi,pid,pidsq,tclasspop,
     1  icost,misscost)
        integer cl(nsample),tclasspop(nclass)
        real pi(nclass),pid(nclass),pidsq(nclass),
     1  misscost(nclass,nclass)
        
        call zerv(tclasspop,nclass)
        do 10 n=1, nsample
        j=cl(n)
        if(j.lt.1.or.j.gt.nclass) then
        write(*,*) j,nclass,n,cl(n)
        stop
        end if
        tclasspop(j)=tclasspop(j)+1
10      continue

        if (ipi.eq.1) then
        do 20 j=1,nclass
        pi(j)=real(tclasspop(j))/nsample
20      continue
        endif
        
        

        do 30 j=1,nclass
        if(tclasspop(j).ge.1) then
        pid(j)=pi(j)/tclasspop(j)
        else
        pid(j)=0
        end if
        pidsq(j)=pid(j)*pid(j)
30      continue

        if(icost.eq.1) then
                do 40 j=1,nclass
                do 50 i=1, nclass
                if (i.eq.j) misscost(i,j)=0
                if (i.ne.j) misscost(i,j)=1
50              continue
40              continue
        endif
        end
        
     
        subroutine testreebag(x,nsample,mdim,treemap,nodestatus,
     1  nodeclass,nrnodes,ndbigtree,nclass,jtree,uts,upper,mbest,rcat,
     1  cat,maxcat,nodex)
        
        
        real x(mdim,nsample),uts(nsample),upper(nrnodes),
     1  rcat(maxcat,nrnodes)
        
        integer treemap(2,nrnodes),cat(mdim),nodex(nsample),
     1  mbest(nrnodes),nodeclass(nrnodes),nodestatus(nrnodes),
     1  jtree(nsample)
    
        call zerv(jtree,nsample)
        call zerv(nodex,nsample)
        
        do n=1,nsample
        kt=1
        do k=1,ndbigtree
        if(nodestatus(kt).eq.-1) then
        jtree(n)=nodeclass(kt)
        nodex(n)=kt
        goto 100
        end if
        mv=mbest(kt)
        if(cat(mv).eq.1) then
        uts(n)=x(mv,n)
        else
        uts(n)=rcat(nint(x(mv,n)),kt)
        end if
                if (uts(n).le.upper(kt)) then 
                        kt=treemap(1,kt)
                else
                        kt=treemap(2,kt)
                endif
        end do
100     continue
        end do
        
        end
        
        
        
        



        subroutine oob(nsample,nclass,jin,cl,jtr,jerr,ncounttr,nout,
     1  errtr,errc,jout,rmargin,q,jest)
        integer jin(nsample),cl(nsample),jtr(nsample),
     1  ncounttr(nclass,nsample),nout(nsample),jerr(nsample),
     1  jest(nsample)

        real rmargin(nsample),q(nclass,nsample)
        
        
        noutc=0
        nmissc=0
        do n=1,nsample
        if(jin(n).eq.0) then
        noutc=noutc+1
        if(jtr(n).ne.cl(n)) nmissc=nmissc+1
        end if
        end do
        errc=real(nmissc)/noutc
        
        nmiss=0
        jout=0
        do n=1,nsample
        if(jin(n).eq.0) then
        ncounttr(jtr(n),n)=ncounttr(jtr(n),n)+1
        nout(n)=nout(n)+1
        end if
        end do
        
        call zerv(jerr,nsample)
        jout=0
        nmiss=0
        do n=1,nsample
        if(nout(n).ge.1) then
        jout=jout+1
        smax=0
        smaxtr=0
        do j=1,nclass
        q(j,n)=real(ncounttr(j,n))/nout(n)
        if(j.ne.cl(n)) smax=amax1(q(j,n),smax)
        if (q(j,n).gt.smaxtr) then
        smaxtr=q(j,n)
        jest(n)=j
        end if
        end do
        if(jest(n).ne.cl(n)) then
        nmiss=nmiss+1
        jerr(n)=1
        end if
        pth=q(cl(n),n)
        rmargin(n)=pth-smax
        end if
        end do
        
        errtr=real(nmiss)/jout
        
        
        end
        
        
        subroutine ndnear(nsample,jin,nodex,prox)       
        real prox(nsample,nsample)
        integer nodex(nsample),jin(nsample)
        do n=1,nsample
        if(jin(n).eq.1) then
        do k=1,nsample
        if(jin(k).eq.0) then
        if(nodex(k).eq.nodex(n)) prox(n,k)=prox(n,k)+1.0
        end if
        end do
        end if
        end do !n
        
        end
                
        
        
        subroutine permobmr(mr,x,tp,tx,jin,nsample,mdim)
        dimension x(mdim,nsample), tp(nsample),tx(nsample),jin(nsample)
        kout=0
        call zervr(tp,nsample)
        do n=1,nsample
        if(jin(n).eq.0) then
        kout=kout+1
        tp(kout)=x(mr,n)
        end if
        end do !n
        call perm1(kout,nsample,tp)
        iout=0
        do n=1,nsample
        tx(n)=x(mr,n)
        if(jin(n).eq.0) then
        iout=iout+1
        x(mr,n)=tp(iout)
        end if
        end do
        end
        
                
        subroutine contimp(mr,x,jin,jts,tx,ncountimp,nclass,nsample,
     1  mdim)
        dimension x(mdim,nsample),tx(nsample),jin(nsample),
     1  ncountimp(nclass,nsample,mdim),jts(nsample)
        do n=1,nsample
        x(mr,n)=tx(n)
        end do
        do n=1,nsample
        if(jin(n).eq.0) then
        ncountimp(jts(n),n,mr)=ncountimp(jts(n),n,mr)+1
        end if
        end do
        end
        
        subroutine finishimp(missimp,ncountimp,nout,cl,nclass,mdim,
     1  nsample,errimp,jout)
        integer missimp(mdim),ncountimp(nclass,nsample,mdim),
     1  nout(nsample),cl(nsample)
        real errimp(mdim) 
        call zerv(missimp,mdim)
        call zervr(errimp,mdim)
        do m1=1,mdim
        do n=1,nsample
        if(nout(n).ge.1) then
        smax=0
        do j=1,nclass
        if (real(ncountimp(j,n,m1))/nout(n).gt.smax) then
        smax=real(ncountimp(j,n,m1))/nout(n)
        imax=j
        end if  
        end do !j
        if(imax.ne.cl(n)) missimp(m1)=missimp(m1)+1
        end if
        end do !n
        end do
        do m1=1,mdim
        errimp(m1)=real(missimp(m1))/jout
        end do
        end
        
	subroutine comptserr(ncount,jtts,jmaxts,clts,ntest,nclass,errts,ltest)
	integer ncount(nclass,ntest),jtts(ntest),clts(ntest),jmaxts(ntest)
	
	missts=0
     	do n=1,ntest
	ncount(jtts(n),n)=ncount(jtts(n),n)+1 
	end do
	do n=1,ntest
	cmax=0
	do j=1,nclass
	if (ncount(j,n).gt.cmax) then
	jmaxts(n)=j
	cmax=ncount(j,n)
	end if
	end do
	end do
	if (ltest.eq.2) then
	nmiss=0
	do n=1,ntest
	if(jmaxts(n).ne.clts(n)) missts=missts+1
	end do
	errts=real(missts)/ntest
	end if
	end
	
      
        
c       MISCELLANOUS SMALL SUBROUTINES  

        
        subroutine zerv(ix,m1)
        integer ix(m1)
        do 10 n=1,m1
        ix(n)=0
10      continue
        end
        
        subroutine zervr(rx,m1)
        real rx(m1)
        do 10 n=1,m1
        rx(n)=0
10      continue
        end
        
        subroutine zerm(mx,m1,m2)
        integer mx(m1,m2)
        do 10 i=1,m1
        do 20 j=1,m2
        mx(i,j)=0
20      continue
10      continue
        end
        
        subroutine zermr(rx,m1,m2)
        real  rx(m1,m2)
        do 10 i=1,m1
        do 20 j=1,m2
        rx(i,j)=0
20      continue
10      continue
        end
        
        
        
        
        
        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
      dimension iperm(kk),v(kk),iu(32),il(32)
      integer t,tt
     
      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
      
        integer function nselect(ns,sm)
        real sm(ns)
        u=randn(1)
        kp=ns
        km=0
        do j=1,10000
        kt=nint(real(kp+km)/2)
        if (u.lt.sm(kt)) then
        kp=kt
        else
        km=kt
        end if
        if (kp-km.eq.1) then
        nselect=kp
	goto 69
        end if
        end do
69      return
        end
        

        
        subroutine perm(ns,ntp)
        integer ntp(ns)
        do 1 n= 1,ns
        ntp(n)=n
1       continue        
        j=ns
11      rnd = randn(1)
        k=int(j*rnd)+1
        jx=ntp(j)
        ntp(j)=ntp(k)
        ntp(k)=jx
        j=j-1
        if(j.gt.1) go to 11
        end
        
        subroutine perm1(np,ns,tp)
        real tp(ns)
        j=np
11      rnd = randn(1)
        k=int(j*rnd)+1
        tx=tp(j)
 210    tp(j)=tp(k)
        tp(k)=tx
        j=j-1
        if(j.gt.1) go to 11
        end

        
        subroutine catran(lc,jc,maxcat)
        dimension lc(maxcat)
        k=int(randn(1)*(jc-1))+1
        do i=1,k
        lc(i)=1
        end do
        do i=k+1,jc
        lc(i)=0
        end do
        j=jc
11      rnd = randn(1)
        k=int(j*rnd)+1
        jx=lc(j)
        lc(j)=lc(k)
        lc(k)=jx
        j=j-1
        if(j.gt.1) go to 11
        end
        
        real function randn(j)
        double precision dseed
        real u
        save dseed
        data dseed /17395/
        call lrnd(dseed,u)
        randn=u
        end

        subroutine lrnd(dseed,u)
        real u
        double precision dseed, d31m1
        data d31m1 /2147483647/
        dseed=dmod(16087*dseed,d31m1)
        u=dseed/d31m1
        return
        end