classtree.c

是基于linux系统的C++程序

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

#include <R.h>
#include <Rmath.h>
#include "rf.h"

#ifdef C_CLASSTREE
void classTree(int *a, int *b, int *class, int *cat, int mdim, int nsample, 
               int nclass, int *treemap, int *bestvar, double *bestsplit,
               double *bestsplitnext, double *tgini, int *nodeStatus,
               int *nodePop, int *nodeStart, double *tclassPop, int maxNodes,
               int nodeSize, int *ncase, int *inBag, int mTry, int *varUsed,
               int *nodeClass, int *treeSize, double *win) {
/* Buildtree consists of repeated calls to two subroutines, Findbestsplit
   and Movedata.  Findbestsplit does just that--it finds the best split of
   the current node.  Movedata moves the data in the split node right and
   left so that the data corresponding to each child node is contiguous.
   The buildtree bookkeeping is different from that in Friedman's original
   CART program.  ncur is the total number of nodes to date.
   nodeStatus(k)=1 if the kth node has been split.  nodeStatus(k)=2 if the
   node exists but has not yet been split, and =-1 of the node is terminal.
   A node is terminal if its size is below a threshold value, or if it is
   all one class, or if all the x-values are equal.  If the current node k
   is split, then its children are numbered ncur+1 (left), and
   ncur+2(right), ncur increases to ncur+2 and the next node to be split is
   numbered k+1.  When no more nodes can be split, buildtree returns to the
   main program.
*/
/*
  integer a(mdim,nsample),cl(nsample),cat(mdim),
  treemap(2,numNodes),bestvar(numNodes),
          bestsplit(numNodes), nodeStatus(numNodes),ta(nsample),
          nodePop(numNodes),nodeStart(numNodes),
          bestsplitnext(numNodes),idmove(nsample),
          ncase(nsample),parent(numNodes),b(mdim,nsample),
          jin(nsample),iv(mred),nodeclass(numNodes),mind(mred)
      
      
      double precision tclasspop(nclass),classpop(nclass,numNodes),
     1     tclasscat(nclass,32),win(nsample),wr(nclass),wc(nclass),
     1     wl(nclass),tgini(mdim), xrand
 */          
    int msplit = 0, i, j;      
    zeroInt(nodeStatus, maxNodes);
    zeroInt(nodeStart, maxNodes);
    zeroInt(nodePop, maxNodes);
    zeroDouble(classPop, nclass * maxNodes);
      
    for (i = 0; i < nclass; ++i) classPop[i] = tclassPop[i];
    ncur = 1;
    nodeStart[0] = 1;
    nodePop[0] = *nuse;
    nodeStatus[0] = NODE_TOSPLIT;
    /* 2: not split yet, 1: split, -1: terminal */      
    /* start main loop */
    for (i = 0; i < numNodes; ++i) {
        if (i > ncur - 1) break;
        if (nodeStatus[i] != NODE_TOSPLIT) continue; 
        /* initialize for next call to findbestsplit */
        ndstart = nodeStart[i];
        ndend = ndstart + nodePop[i] - 1;
        for (j = 0; j < nclass; ++j) {
            tclassPop[j] = classPop[j + i * nclass];
        }
        jstat = 0;
        F77_CALL(findbestsplit)(a, b, cl, mdim, nsample, nclass, cat,
                                ndstart, ndend, tclassPop, tclasscat,
                                &msplit, &decsplit, &nbest, ncase, &jstat,
                                inBag, mTry, win, wr, wc, wl, mred, i, mind);
        if (jstat == 1) { 
            nodeStatus[i] = NODE_TERMINAL;
            continue;
        } else { 
            bestvar[i] = msplit;
            varUsed[msplit - 1] = 1;
            tgini[msplit - 1] += decsplit;
            if (cat[msplit-1] == 1) {
                bestsplit[i] = a[msplit - 1  + nbest * mdim];
                bestsplitnext[i] = a[msplit - 1 + (nbest + 1) * mdim];
            } else {
                  bestsplit[i] = nbest;
                  bestsplitnext[i] = 0;
            }
        }          
        F77_CALL(movedata)(a, ta, mdim, nsample, ndstart, ndend, idmove,
                           ncase, msplit, cat, nbest, ndendl);
        /* 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;
        /* find class populations in both nodes */
        for (n = ndstart; n <= ndendl; ++n) {
            nc = ncase[n];
            j = class[nc-1];
            classPop[j - 1 + (ncur+1)*mdim] += win[nc - 1];
        }
        for (n = ndendl + 1; n <= ndend; ++n) {
            nc = ncase[n];
            j = cl[nc - 1];
            classPop[j - 1 + (ncur+2) * mdim] += win[nc - 1];
        }
        /* check on nodeStatus */
        nodeStatus[ncur + 1] = NODE_TOSPLIT;
        nodeStatus[ncur + 2] = NODE_TOSPLIT;
        if (nodePop[ncur + 1] <= ndsize) nodeStatus[ncur+1] = NODE_TERMINAL;
        if (nodePop[ncur + 2] <= ndsize) nodeStatus[ncur+2] = NODE_TERMINAL;
        popt1 = 0;
        popt2 = 0;
        for (j = 0; j < nclass; ++j) {
            popt1 += classPop[j + (ncur+1) * mdim];
            popt2 += classPop[j + (ncur+2) * mdim];
        }
        for (j = 0; j < nclass; ++j) {
            if (classPop[j + (ncur+1) * mdim] == popt1) 
                nodeStatus[ncur+1] = NODE_TERMINAL;
            if (classPop[j + (ncur+2) * mdim] == popt2) 
                nodeStatus[ncur+2] = NODE_TERMINAL;
        }
        
        treemap[i * 2] = ncur + 1;
        treemap[1 + i * 2] = ncur + 2;
        nodeStatus[i] = NODE_INTERIOR;
        ncur += 2;
        if (ncur >= numNodes) break;
    }
    ndbigtree = numNodes;
    for (k = numNodes-1; k >= 0; --k) {
        if (nodeStatus[k] == 0) ndbigtree--;
        if (nodeStatus[k] == NODE_TOSPLIT) nodeStatus[k] = NODE_TERMINAL;
    }
    for (k = 0; k < ndbigtree; ++k) {
        if (nodeStatus[k] == NODE_TERMINAL) {
            pp = 0;
            for (j = 0; j < nclass; ++j) {
                if (classPop[j + k * nclass] > pp) {
                    nodeClass[k] = j;
                    pp = classPop[j + k * nclass];
                }
                /* Break ties at random: */
                if (classPop[j + k * nclass] == pp && unif_rand() > 0.5) {
                    nodeClass[k] = j;
                    pp = classPop[j + k * nclass];
                }
            }
        }
    }
}



void findBestSplit(int *a, double *b, int *class, int mDim, int nSample, 
                   int nClass, int *nCat, int maxCat, int ndStart, int ndEnd, 
                   double *classCount, double *classCatTable, 
                   int *splitVar, double *decGini, int *bestSplit,
                   int *ncase, int *splitStatus, int *inBag, int mtry,
                   double *weight, double *wr, double *wc, double *wl,
                   int *currentNode, int *mind) {
/*
      subroutine findbestsplit(a, b, cl, mdim, nsample, nclass, cat,
     1     maxcat, ndstart, ndend, tclasspop, tclasscat, msplit, 
     2     decsplit, nbest, ncase, jstat, jin, mtry, win, wr, wc, wl,
     3     mred, kbuild, mind) */
/*
     For the best split, msplit is the variable split on. decsplit is the
     dec. in impurity.  If msplit is numerical, nsplit is the case number
     of value of msplit split on, and nsplitnext is the case number of the
     next larger value of msplit.  If msplit is categorical, then nsplit is
     the coding into an integer of the categories going left.
*/


    integer a(mdim,nsample), cl(nsample), cat(mdim),
     1     ncase(nsample), b(mdim,nsample), jin(nsample), nn, j          
      double precision tclasspop(nclass), tclasscat(nclass,32), dn(32),
     1     win(nsample), wr(nclass), wc(nclass), wl(nclass), xrand
      integer mind(mred), ncmax, ncsplit,nhit
        ncmax = 10;
    ncsplit = 512;
    /* compute initial values of numerator and denominator of Gini */
    parentNum = 0.0;
    parentDen = 0.0;
    for (i = 0; i < nClass; ++i) {
        parentNum += classCount[i] * classCount[i];
        parentDen += classCount[i];
    }
    crit0 = pno / pdo;
    *splitStatus = 0;
    critmax = -1.0e25;
    for (i = 0; i < mDim; ++i) mind[i] = i;

    /* start main loop through variables to find best split. */
    last = mDim - 1;
    for (i = 0, i < mtry; ++i) {
        /* sample mtry variables w/o replacement. */
        j = (int) (unif_rand() * (last + 1));
        mvar = mIndex[j];
        swapInt(mIndex[j], mIndex[last]);
        last--;

        lcat = nCat[mvar];
        if (lcat == 1) {
            /* Split on a numerical predictor. */
            rightNum = parentNum;
            rightDen = parentDen;
            leftNum = 0.0;
            leftDen = 0.0;
            zeroDouble(wl, nClass);
            for (j = 0; j < nClass; ++j) wr[j] = classCount[j];
	    ntie = 1;
            for (j = ndstart; j <= ndend - 1; ++j) {
                nc = a[mvar, j-1];
                u = weight[nc];
                k = class[nc];