Tree.cc

Go to the documentation of this file.

//                            Tree.cxx                                  //
// =====================================================================//
// This is the object implementation of a decision tree.                //
// References include                                                   //
//    *Elements of Statistical Learning by Hastie,                      //
//     Tibshirani, and Friedman.                                        //
//    *Greedy Function Approximation: A Gradient Boosting Machine.      //
//     Friedman. The Annals of Statistics, Vol. 29, No. 5. Oct 2001.    //
//    *Inductive Learning of Tree-based Regression Models. Luis Torgo.  //
//                                                                      //

// _______________________Includes_______________________________________//

#include "L1Trigger/L1TMuonEndCap/interface/bdt/Tree.h"

#include <iostream>
#include <sstream>
#include <cmath>

// _______________________Constructor(s)________________________________//

using namespace emtf;

Tree::Tree()
{
    rootNode = new Node("root");

    terminalNodes.push_back(rootNode);
    numTerminalNodes = 1;
    boostWeight = 0;
    xmlVersion = 2017;
}

Tree::Tree(std::vector< std::vector<Event*> >& cEvents)
{
    rootNode = new Node("root");
    rootNode->setEvents(cEvents);

    terminalNodes.push_back(rootNode);
    numTerminalNodes = 1;
    boostWeight = 0;
    xmlVersion = 2017;
}
// _______________________Destructor____________________________________//


Tree::~Tree()
{
// When the tree is destroyed it will delete all of the nodes in the tree.
// The deletion begins with the rootnode and continues recursively.
   if(rootNode) delete rootNode;
}

Tree::Tree(const Tree &tree)
{
    // unfortunately, authors of these classes didn't use const qualifiers
    rootNode = copyFrom(const_cast<Tree&>(tree).getRootNode());
    numTerminalNodes = tree.numTerminalNodes;
    rmsError = tree.rmsError;
    boostWeight = tree.boostWeight;
    xmlVersion = tree.xmlVersion;

    terminalNodes.resize(0);
    // find new leafs
    findLeafs(rootNode,terminalNodes);

}

Tree& Tree::operator=(const Tree &tree){
    if(rootNode) delete rootNode;
    // unfortunately, authors of these classes didn't use const qualifiers
    rootNode = copyFrom(const_cast<Tree&>(tree).getRootNode());
    numTerminalNodes = tree.numTerminalNodes;
    rmsError = tree.rmsError;
    boostWeight = tree.boostWeight;
    xmlVersion = tree.xmlVersion;

    terminalNodes.resize(0);
    // find new leafs
    findLeafs(rootNode,terminalNodes);


    return *this;
}

Node* Tree::copyFrom(const Node *local_root)
{
    // end-case
    if( !local_root ) return nullptr;

    Node *lr = const_cast<Node*>(local_root);

    // recursion
    Node *left_new_child  = copyFrom( lr->getLeftDaughter()  );
    Node *right_new_child = copyFrom( lr->getRightDaughter() );

    // performing main work at this level
    Node *new_local_root  = new Node( lr->getName() );
    if( left_new_child  ) left_new_child ->setParent(new_local_root);
    if( right_new_child ) right_new_child->setParent(new_local_root);
    new_local_root->setLeftDaughter ( left_new_child  );
    new_local_root->setRightDaughter( right_new_child );
    new_local_root->setErrorReduction( lr->getErrorReduction() );
    new_local_root->setSplitValue( lr->getSplitValue() );
    new_local_root->setSplitVariable( lr->getSplitVariable() );
    new_local_root->setFitValue( lr->getFitValue() );
    new_local_root->setTotalError( lr->getTotalError() );
    new_local_root->setAvgError( lr->getAvgError() );
    new_local_root->setNumEvents( lr->getNumEvents() );
//    new_local_root->setEvents( lr->getEvents() ); // no ownership assumed for the events anyways

    return new_local_root;
}

void Tree::findLeafs(Node *local_root, std::list<Node*> &tn)
{
    if( !local_root->getLeftDaughter() && !local_root->getRightDaughter() ){
        // leaf or ternimal node found
        tn.push_back(local_root);
        return;
    }

    if( local_root->getLeftDaughter() )
        findLeafs( local_root->getLeftDaughter(), tn );

    if( local_root->getRightDaughter() )
        findLeafs( local_root->getRightDaughter(), tn );
}

Tree::Tree(Tree && tree)
{
  if(rootNode) delete rootNode; // this line is the only reason not to use default move constructor
  rootNode = tree.rootNode;
  terminalNodes = std::move(tree.terminalNodes);
  numTerminalNodes = tree.numTerminalNodes;
  rmsError = tree.rmsError;
  boostWeight = tree.boostWeight;
  xmlVersion = tree.xmlVersion;
}

// ______________________Get/Set________________________________________//

void Tree::setRootNode(Node *sRootNode)
{
    rootNode = sRootNode;
}

Node * Tree::getRootNode()
{
     return rootNode;
}

// ----------------------------------------------------------------------

void Tree::setTerminalNodes(std::list<Node*>& sTNodes)
{
    terminalNodes = sTNodes;
}

std::list<Node*>& Tree::getTerminalNodes()
{
    return terminalNodes;
}

// ----------------------------------------------------------------------

int Tree::getNumTerminalNodes()
{
    return numTerminalNodes;
}

// ______________________Performace_____________________________________//

void Tree::calcError()
{
// Loop through the separate predictive regions (terminal nodes) and
// add up the errors to get the error of the entire space.

    double totalSquaredError = 0;

    for(std::list<Node*>::iterator it=terminalNodes.begin(); it!=terminalNodes.end(); it++)
    {
        totalSquaredError += (*it)->getTotalError();
    }
    rmsError = sqrt( totalSquaredError/rootNode->getNumEvents() );
}

// ----------------------------------------------------------------------

void Tree::buildTree(int nodeLimit)
{
    // We greedily pick the best terminal node to split.
    double bestNodeErrorReduction = -1;
    Node* nodeToSplit = nullptr;

    if(numTerminalNodes == 1)
    {
        rootNode->calcOptimumSplit();
        calcError();
//        std::cout << std::endl << "  " << numTerminalNodes << " Nodes : " << rmsError << std::endl;
    }

    for(std::list<Node*>::iterator it=terminalNodes.begin(); it!=terminalNodes.end(); it++)
    {
       if( (*it)->getErrorReduction() > bestNodeErrorReduction )
       {
           bestNodeErrorReduction = (*it)->getErrorReduction();
           nodeToSplit = (*it);
       }
    }

    //std::cout << "nodeToSplit size = " << nodeToSplit->getNumEvents() << std::endl;

    // If all of the nodes have one event we can't add any more nodes and reduce the error.
    if(nodeToSplit == nullptr) return;

    // Create daughter nodes, and link the nodes together appropriately.
    nodeToSplit->theMiracleOfChildBirth();

    // Get left and right daughters for reference.
    Node* left = nodeToSplit->getLeftDaughter();
    Node* right = nodeToSplit->getRightDaughter();

    // Update the list of terminal nodes.
    terminalNodes.remove(nodeToSplit);
    terminalNodes.push_back(left);
    terminalNodes.push_back(right);
    numTerminalNodes++;

    // Filter the events from the parent into the daughters.
    nodeToSplit->filterEventsToDaughters();

    // Calculate the best splits for the new nodes.
    left->calcOptimumSplit();
    right->calcOptimumSplit();

    // See if the error reduces as we add more nodes.
    calcError();

    if(numTerminalNodes % 1 == 0)
    {
//        std::cout << "  " << numTerminalNodes << " Nodes : " << rmsError << std::endl;
    }

    // Repeat until done.
    if(numTerminalNodes <  nodeLimit) buildTree(nodeLimit);
}

// ----------------------------------------------------------------------

void Tree::filterEvents(std::vector<Event*>& tEvents)
{
// Use trees which have already been built to fit a bunch of events
// given by the tEvents vector.

    // Set the events to be filtered.
    rootNode->getEvents() = std::vector< std::vector<Event*> >(1);
    rootNode->getEvents()[0] = tEvents;

    // The tree now knows about the events it needs to fit.
    // Filter them into a predictive region (terminal node).
    filterEventsRecursive(rootNode);
}

// ----------------------------------------------------------------------

void Tree::filterEventsRecursive(Node* node)
{
// Filter the events repeatedly into the daughter nodes until they
// fall into a terminal node.

    Node* left = node->getLeftDaughter();
    Node* right = node->getRightDaughter();

    if(left == nullptr || right == nullptr) return;

    node->filterEventsToDaughters();

    filterEventsRecursive(left);
    filterEventsRecursive(right);
}

// ----------------------------------------------------------------------

Node* Tree::filterEvent(Event* e)
{
// Use trees which have already been built to fit a bunch of events
// given by the tEvents vector.

    // Filter the event into a predictive region (terminal node).
    Node* node = filterEventRecursive(rootNode, e);
    return node;
}

// ----------------------------------------------------------------------

Node* Tree::filterEventRecursive(Node* node, Event* e)
{
// Filter the event repeatedly into the daughter nodes until it
// falls into a terminal node.


    Node* nextNode = node->filterEventToDaughter(e);
    if(nextNode == nullptr) return node;

    return filterEventRecursive(nextNode, e);
}

// ----------------------------------------------------------------------


void Tree::rankVariablesRecursive(Node* node, std::vector<double>& v)
{
// We recursively go through all of the nodes in the tree and find the
// total error reduction for each variable. The one with the most
// error reduction should be the most important.

    Node* left = node->getLeftDaughter();
    Node* right = node->getRightDaughter();

    // Terminal nodes don't contribute to error reduction.
    if(left==nullptr || right==nullptr) return;

    int sv =  node->getSplitVariable();
    double er = node->getErrorReduction();

    //if(sv == -1)
    //{
      //std::cout << "ERROR: negative split variable for nonterminal node." << std::endl;
      //std::cout << "rankVarRecursive Split Variable = " << sv << std::endl;
      //std::cout << "rankVarRecursive Error Reduction = " << er << std::endl;
    //}

    // Add error reduction to the current total for the appropriate
    // variable.
    v[sv] += er;

    rankVariablesRecursive(left, v);
    rankVariablesRecursive(right, v);

}

// ----------------------------------------------------------------------

void Tree::rankVariables(std::vector<double>& v)
{
    rankVariablesRecursive(rootNode, v);
}

// ----------------------------------------------------------------------


void Tree::getSplitValuesRecursive(Node* node, std::vector<std::vector<double>>& v)
{
// We recursively go through all of the nodes in the tree and find the
// split points used for each split variable.

    Node* left = node->getLeftDaughter();
    Node* right = node->getRightDaughter();

    // Terminal nodes don't contribute.
    if(left==nullptr || right==nullptr) return;

    int sv =  node->getSplitVariable();
    double sp = node->getSplitValue();

    if(sv == -1)
    {
        std::cout << "ERROR: negative split variable for nonterminal node." << std::endl;
        std::cout << "rankVarRecursive Split Variable = " << sv << std::endl;
    }

    // Add the split point to the list for the correct split variable.
    v[sv].push_back(sp);

    getSplitValuesRecursive(left, v);
    getSplitValuesRecursive(right, v);

}

// ----------------------------------------------------------------------

void Tree::getSplitValues(std::vector<std::vector<double>>& v)
{
    getSplitValuesRecursive(rootNode, v);
}

// ______________________Storage/Retrieval______________________________//

template <typename T>
std::string numToStr( T num )
{
// Convert a number to a string.
    std::stringstream ss;
    ss << num;
    std::string s = ss.str();
    return  s;
}

// ----------------------------------------------------------------------

void Tree::addXMLAttributes(TXMLEngine* xml, Node* node, XMLNodePointer_t np)
{
    // Convert Node members into XML attributes
    // and add them to the XMLEngine.
    xml->NewAttr(np, nullptr, "splitVar", numToStr(node->getSplitVariable()).c_str());
    xml->NewAttr(np, nullptr, "splitVal", numToStr(node->getSplitValue()).c_str());
    xml->NewAttr(np, nullptr, "fitVal", numToStr(node->getFitValue()).c_str());
}

// ----------------------------------------------------------------------

void Tree::saveToXML(const char* c)
{

    TXMLEngine* xml = new TXMLEngine();

    // Add the root node.
    XMLNodePointer_t root = xml->NewChild(nullptr, nullptr, rootNode->getName().c_str());
    addXMLAttributes(xml, rootNode, root);

    // Recursively write the tree to XML.
    saveToXMLRecursive(xml, rootNode, root);

    // Make the XML Document.
    XMLDocPointer_t xmldoc = xml->NewDoc();
    xml->DocSetRootElement(xmldoc, root);

    // Save to file.
    xml->SaveDoc(xmldoc, c);

    // Clean up.
    xml->FreeDoc(xmldoc);
    delete xml;
}

// ----------------------------------------------------------------------

void Tree::saveToXMLRecursive(TXMLEngine* xml, Node* node, XMLNodePointer_t np)
{
    Node* l = node->getLeftDaughter();
    Node* r = node->getRightDaughter();

    if(l==nullptr || r==nullptr) return;

    // Add children to the XMLEngine.
    XMLNodePointer_t left = xml->NewChild(np, nullptr, "left");
    XMLNodePointer_t right = xml->NewChild(np, nullptr, "right");

    // Add attributes to the children.
    addXMLAttributes(xml, l, left);
    addXMLAttributes(xml, r, right);

    // Recurse.
    saveToXMLRecursive(xml, l, left);
    saveToXMLRecursive(xml, r, right);
}

// ----------------------------------------------------------------------

void Tree::loadFromXML(const char* filename)
{
    // First create the engine.
    TXMLEngine* xml = new TXMLEngine;

    // Now try to parse xml file.
    XMLDocPointer_t xmldoc = xml->ParseFile(filename);
    if (xmldoc==nullptr)
    {
        delete xml;
        return;
    }

    // Get access to main node of the xml file.
    XMLNodePointer_t mainnode = xml->DocGetRootElement(xmldoc);

    // the original 2016 pT xmls define the source tree node to be the top-level xml node
    // while in 2017 TMVA's xmls every decision tree is wrapped in an extra block specifying boostWeight parameter
    // I choose to identify the format by checking the top xml node name that is a "BinaryTree" in 2017
    if( std::string("BinaryTree") == xml->GetNodeName(mainnode) ){
        XMLAttrPointer_t attr = xml->GetFirstAttr(mainnode);
        attr = xml->GetNextAttr(attr);
        boostWeight = (attr ? strtod(xml->GetAttrValue(attr),nullptr) : 0);
        // step inside the top-level xml node
        mainnode = xml->GetChild(mainnode);
        xmlVersion = 2017;
    } else {
        boostWeight = 0;
        xmlVersion = 2016;
    }
    // Recursively connect nodes together.
    loadFromXMLRecursive(xml, mainnode, rootNode);

    // Release memory before exit
    xml->FreeDoc(xmldoc);
    delete xml;
}

// ----------------------------------------------------------------------

void Tree::loadFromXMLRecursive(TXMLEngine* xml, XMLNodePointer_t xnode, Node* tnode)
{

    // Get the split information from xml.
    XMLAttrPointer_t attr = xml->GetFirstAttr(xnode);
    std::vector<std::string> splitInfo(3);
    if( xmlVersion >= 2017 ){
        for(unsigned int i=0,j=0; i<10; i++)
        {
          if(i==3 || i==4 || i==6){
              splitInfo[j++] = xml->GetAttrValue(attr);
          }
          attr = xml->GetNextAttr(attr);
        }
    } else {
        for(unsigned int i=0; i<3; i++)
        {
            splitInfo[i] = xml->GetAttrValue(attr);
            attr = xml->GetNextAttr(attr);
        }
    }

    // Convert strings into numbers.
    std::stringstream converter;
    int splitVar;
    double splitVal;
    double fitVal;

    converter << splitInfo[0];
    converter >> splitVar;
    converter.str("");
    converter.clear();

    converter << splitInfo[1];
    converter >> splitVal;
    converter.str("");
    converter.clear();

    converter << splitInfo[2];
    converter >> fitVal;
    converter.str("");
    converter.clear();

    // Store gathered splitInfo into the node object.
    tnode->setSplitVariable(splitVar);
    tnode->setSplitValue(splitVal);
    tnode->setFitValue(fitVal);

    // Get the xml daughters of the current xml node.
    XMLNodePointer_t xleft = xml->GetChild(xnode);
    XMLNodePointer_t xright = xml->GetNext(xleft);

    // If there are no daughters we are done.
    if(xleft == nullptr || xright == nullptr) return;

    // If there are daughters link the node objects appropriately.
    tnode->theMiracleOfChildBirth();
    Node* tleft = tnode->getLeftDaughter();
    Node* tright = tnode->getRightDaughter();

    // Update the list of terminal nodes.
    terminalNodes.remove(tnode);
    terminalNodes.push_back(tleft);
    terminalNodes.push_back(tright);
    numTerminalNodes++;

    loadFromXMLRecursive(xml, xleft, tleft);
    loadFromXMLRecursive(xml, xright, tright);
}

void Tree::loadFromCondPayload(const L1TMuonEndCapForest::DTree& tree)
{
    // start fresh in case this is not the only call to construct a tree
    if( rootNode ) delete rootNode;
    rootNode = new Node("root");

    const L1TMuonEndCapForest::DTreeNode& mainnode = tree[0];
    loadFromCondPayloadRecursive(tree, mainnode, rootNode);
}

void Tree::loadFromCondPayloadRecursive(const L1TMuonEndCapForest::DTree& tree, const L1TMuonEndCapForest::DTreeNode& node, Node* tnode)
{
    // Store gathered splitInfo into the node object.
    tnode->setSplitVariable(node.splitVar);
    tnode->setSplitValue(node.splitVal);
    tnode->setFitValue(node.fitVal);

    // If there are no daughters we are done.
    if( node.ileft == 0 || node.iright == 0) return; // root cannot be anyone's child
    if( node.ileft  >= tree.size() ||
        node.iright >= tree.size() ) return; // out of range addressing on purpose

    // If there are daughters link the node objects appropriately.
    tnode->theMiracleOfChildBirth();
    Node* tleft = tnode->getLeftDaughter();
    Node* tright = tnode->getRightDaughter();

    // Update the list of terminal nodes.
    terminalNodes.remove(tnode);
    terminalNodes.push_back(tleft);
    terminalNodes.push_back(tright);
    numTerminalNodes++;

    loadFromCondPayloadRecursive(tree, tree[node.ileft], tleft);
    loadFromCondPayloadRecursive(tree, tree[node.iright], tright);
}