probabilisticregressiontreenode.h

Go to the documentation of this file.

/*

Copyright (c) 2003, Cornell University
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

   - Redistributions of source code must retain the above copyright notice,
       this list of conditions and the following disclaimer.
   - Redistributions in binary form must reproduce the above copyright
       notice, this list of conditions and the following disclaimer in the
       documentation and/or other materials provided with the distribution.
   - Neither the name of Cornell University nor the names of its
       contributors may be used to endorse or promote products derived from
       this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.

*/

// -*- C++ -*-

#if !defined _CLUS_PROBABILISTICREGRESSIONTREENODE_H_
#define _CLUS_PROBABILISTICREGRESSIONTREENODE_H_

#include "vec.h"
#include "distribution.h"
#include "dynamicbuffer.h"

using namespace TNT;

#define EMMAXTRIALS 3

namespace CLUS
{

/** Class used in building regression trees.
    Every class is a node in the tree, there can be 0 or 2 children
    NodeId determines the position of the node in the tree. The nodeId
    od children is alwais 2*nodeId and 2*nodeId+1. The nodeId starts from 1
    for the root. 
 
    T_Distribution is the distribution that approximates the data
    T_Regressor is the regressor that commes with T_Distribution
    T_Spliter is the class that in a final scan can compute the 
    splitting predicate
*/
template< class T_Distribution, class T_Regressor, class T_Splitter >
class BinaryProbabilisticRegressionTreeNode
{
protected:
    /// unique identifier of the cluster for a regression tree.
    /// the id of the children is always 2*Id, 2*Id+1
    int nodeId;

    /// how many times EM was attempted
    int emTrials;

    /// the number of continuous split variables
    int csDim;

    /// the number of regressors
    int regDim;
    
    /// the state of the node. At creation em. At load stable
    enum state { stable, em, split, regression } State;
    
    /// the children of this node
    BinaryProbabilisticRegressionTreeNode< T_Distribution, T_Regressor, T_Splitter > * Children[2];
    
    /// splitter for split criterion
    T_Splitter Splitter;
    
    /// regressor for the node
    T_Regressor Regressor;
    
    /// distributions for the children for the learning process    
    T_Distribution* Distributions[2];
    
    /// the distribution of the parent. Usefull in the learning process to pick initial values for the new distributions
    T_Distribution* parentDistribution;

    /// keep the EM samples here temporarily
    DynamicBuffer* buffer; 

    /// sum of squared differences between prediction and true value
    double pruningCost;

    /// number of samples in this node for pruning
    int pruningSamples; 

    void RandomDistributions(void)
    {
        /*      if (parentDistribution != 0){
         Distributions[0]->RandomDistribution(2,*parentDistribution);
         Distributions[1]->RandomDistribution(2,*parentDistribution);
         } else { 
        */

        // Distributions are build on normalized data
        Distributions[0]->RandomDistribution(2);
        Distributions[1]->RandomDistribution(2);

        // }
    }

    /**
    Perform an EM step on the data in buff and return the convergence factor 
    and compute the likelihood. 
    */
    double EMStep(double& Likelihood)
    {
        Likelihood = 0.0;
        for (double* X=buffer->begin(); X<buffer->end(); X+=regDim+2)
        {
            // Do the EM learning on the two distributions
            double probability=X[regDim+1];

            double p0=Distributions[0]->LearnProbability(X);
            double p1=Distributions[1]->LearnProbability(X);
            Likelihood+=log((p0+p1)/probability);
            Distributions[0]->NormalizeLearnProbability((p0+p1)/probability,2);
            Distributions[1]->NormalizeLearnProbability((p0+p1)/probability,2);
        }

        double c0=Distributions[0]->UpdateParameters();
        double c1=Distributions[1]->UpdateParameters();

        return ( (c0+c1)/2.0 );
    }


public:
    /// Constructor for a leaf
    BinaryProbabilisticRegressionTreeNode(int NodeId, int CsDim, T_Regressor& regressor,  T_Distribution* ParentDistribution = 0):
            nodeId(NodeId), csDim(CsDim), Splitter(), Regressor(regressor), parentDistribution(ParentDistribution)
    {
        Children[0]=Children[1]=0;
        Distributions[0]=Distributions[0]=0;
        State=stable; // do no learning on this node
        emTrials=0;
        buffer=0;
    }

    /// Constructor for an intermediate node
    BinaryProbabilisticRegressionTreeNode( int NodeId,
                                           const Vector<int> & DDomainSize,
                                           int CsplitDim, int RegDim,
                                           T_Regressor& regressor,
                                           T_Distribution* ParentDistribution = 0):
            nodeId(NodeId),csDim(CsplitDim), regDim(RegDim), State(em),
            Splitter(DDomainSize,CsplitDim,RegDim), Regressor(regressor),
            parentDistribution(ParentDistribution)
    {

        Children[0]=Children[1]=0;
        Distributions[0] = new T_Distribution(regDim);
        Distributions[1] = new T_Distribution(regDim);
        RandomDistributions();
        emTrials=0;
        buffer=0;
    }

    ~BinaryProbabilisticRegressionTreeNode(void)
    {
        if (Children[0]!=0)
            delete Children[0];
        Children[0]=0;

        if (Children[1]!=0)
            delete Children[1];
        Children[1]=0;

        /*      if (Distributions[0]!=0)
         delete Distributions[0];

         if (Distributions[1]!=0)
         delete Distributions[1]; */
    }

    int GetNodeId(void)
    {
        return nodeId;
    }

    void ComputeSizesTree(int& nodes, int& term_nodes)
    {
        nodes++;
        if (Children[0]!=0 && Children[1]!=0)
        {
            Children[0]->ComputeSizesTree(nodes,term_nodes);
            Children[1]->ComputeSizesTree(nodes,term_nodes);
        }
        else
        {
            term_nodes++;
        }
    }

    void StartLearningEpoch(void)
    {
        switch (State)
        {
        case stable:
            if (Children[0]!=0)
            {
                Children[0]->StartLearningEpoch();
                Children[1]->StartLearningEpoch();
            }
            break;
        case em:
            buffer = new DynamicBuffer(regDim+2); // allocate the temporary buffer. we keep weights here too in last position
            break;
        case split:
            // delete parentDistribution; // we dont' need the parent distribution anymore and we get rid of it
            // parentDistribution=0;
            Splitter.InitializeSplitStatistics();
            break;
        case regression:
            // nothing to do
            break;
        }
    }

    void LearnSample(const int* Dvars, const double* Cvars,
                     double probability, double threshold)
    {
        double p0,p1;

        if (probability<threshold)
            return;

        switch (State)
        {
        case stable:
            if (Children[0]!=0)
            {
                double pChild1=Splitter.ProbabilityFirstBranch(Dvars,Cvars);
                Children[0]->LearnSample(Dvars,Cvars,pChild1*probability, threshold);
                Children[1]->LearnSample(Dvars,Cvars,(1.0-pChild1)*probability, threshold);
            }
            break;
        case em:
            // First normalize, otherwise em converges to wrong things (no substructure identified)
            // The normalization makes transforms the data to be centered on 0 and have identity
            // spartity matrix (look superficially like a boal

            /// @todo put some cutoff here to use less data as you go down

            assert(buffer!=0);
            {
                double* cBufferLine=buffer->next();
                parentDistribution->NormalizeData(Cvars+csDim,cBufferLine);
                cBufferLine[regDim+1]=probability;
            }
            break;
        case split:
            p0=Distributions[0]->LearnProbability(Cvars+csDim);
            p1=Distributions[1]->LearnProbability(Cvars+csDim);
            Splitter.UpdateSplitStatistics(Dvars, Cvars, p0/(p0+p1), p1/(p0+p1), probability );
            break;
        case regression:
            {
                double pChild1=Splitter.ProbabilityLeft(Dvars,Cvars);

                if (pChild1 > threshold)
                {
                    p0=Distributions[0]->LearnProbability(Cvars+csDim);
                    Distributions[0]->NormalizeLearnProbability(p0/(pChild1*probability));
                }
                if (1.0-pChild1 > threshold)
                {
                    p1=Distributions[1]->LearnProbability(Cvars+csDim);
                    Distributions[1]->NormalizeLearnProbability(p1/((1.0-pChild1)*probability));
                }
            }
            break;
        }
    }

    /// Returns true if more learning has to be donne in future
    bool StopLearningEpoch(int splitType, int emRestarts, int emMaxIterations,
                           double convergenceLim, int min_no_datapoints)
    {
        bool moresplits=false;
        int emIterations;

        T_Regressor* regressor0=0, * regressor1=0;
        switch (State)
        {
        case stable:
            if (Children[0]!=0)
                return Children[0]->StopLearningEpoch(splitType, emRestarts, emMaxIterations,
                                                      convergenceLim, min_no_datapoints)
                       | Children[1]->StopLearningEpoch(splitType, emRestarts, emMaxIterations,
                                                        convergenceLim, min_no_datapoints); // || cannot be used since is shortcuted
            else
                return false;
        case em:
            emIterations=0;

            {
                cout << "XXXNodeID: " << nodeId << " noDatapoints: " << buffer->dim() << endl;

                // do a number of random restarts and pick the starting point with the maximum likelihood
                T_Distribution best_d0(0);
                T_Distribution best_d1(0);
                double best_Likelihood=-1.0e+100;
                for (int repetition=0; repetition<emRestarts; repetition++)
                {
                    // restart randomly
                    RandomDistributions();
                    double Likelihood;
                    // do two EM steps
                    EMStep(Likelihood);
                    EMStep(Likelihood);

                    cout << "NodeID=" << nodeId << " repetition=" << repetition ;
                    cout << " Likelihood=" << Likelihood << endl;

                    if (Likelihood > best_Likelihood)
                    {
                        best_d0 = *(Distributions[0]);
                        best_d1 = *(Distributions[1]);
                        best_Likelihood = Likelihood;
                    }
                }

                *(Distributions[0])=best_d0;
                *(Distributions[1])=best_d1;
            }

            // continue the learning process with the best starting point

            while (emIterations<emMaxIterations)
            {
                double Likelihood;
                double convFactor = EMStep(Likelihood);

                cout << "NodeID=" << nodeId << " Likelihood=" << Likelihood << endl;

                emIterations++;

                if ( Distributions[0]->HasZeroWeight() || Distributions[1]->HasZeroWeight() )
                {
                    if (emTrials < EMMAXTRIALS)
                    {
                        emTrials++;
                        cerr << "One of the distributions got killed. Starting again" << endl;
                        RandomDistributions();
                        emIterations=0;
                    }
                    else
                    {
                        cerr << "Tried " << EMMAXTRIALS << " times and didn't work. Making the node a leaf." << endl;
                        goto makeleaf;
                    }
                }

                if (!finite(convFactor) || convFactor <= convergenceLim)
                    break;
            }

            // denormalize parameters for the distributions
            Distributions[0]->DenormalizeParameters(parentDistribution);
            Distributions[1]->DenormalizeParameters(parentDistribution);
            State=split;

            delete buffer;
            buffer=0; // to avoid deleting it again

            // print the distribution information for the two distributions
            Distributions[0]->SaveToStream(cout);
            cout << endl;
            Distributions[1]->SaveToStream(cout);
            cout << endl;

            return true;

        case split:

            // cout << "Node: " << nodeId << " being split" << endl;

            State=regression;

            if (Splitter.ComputeSplitVariable(splitType)!=0)
                goto makeleaf;

            Splitter.DeleteTemporaryStatistics();

            return true;
            break;

        case regression:
            cout << "Node:" << nodeId << " finishing regression" << endl;

            Distributions[0]->UpdateParameters();
            Distributions[1]->UpdateParameters();

            regressor0 = dynamic_cast<T_Regressor*> ( Distributions[0]->CreateRegressor() );
            regressor1 = dynamic_cast<T_Regressor*> ( Distributions[1]->CreateRegressor() );

            if ( !regressor0 || !regressor1 )
                goto makeleaf;

            State=stable;

            if (Splitter.MoreSplits(0, min_no_datapoints))
            {
                // more splits in future for child 1, create normal child
                moresplits=true;
                Children[0] = new BinaryProbabilisticRegressionTreeNode< T_Distribution, T_Regressor, T_Splitter >
                              ( nodeId*2, Splitter.GetDDomainSize(), Splitter.GetCSplitDim(),
                                Splitter.GetRegDim(), *regressor0, Distributions[0] );
            }
            else
            {
                // no more splits in future, create leaf children
                Children[0]=new BinaryProbabilisticRegressionTreeNode< T_Distribution, T_Regressor, T_Splitter >
                            ( nodeId*2, csDim, *regressor0, Distributions[0] );
            }

            if (Splitter.MoreSplits(1, min_no_datapoints))
            {
                moresplits=true;
                Children[1] = new BinaryProbabilisticRegressionTreeNode< T_Distribution, T_Regressor, T_Splitter >
                              ( nodeId*2+1, Splitter.GetDDomainSize(), Splitter.GetCSplitDim(),
                                Splitter.GetRegDim(), *regressor1, Distributions[1] );
            }
            else
            {
                Children[1]=new BinaryProbabilisticRegressionTreeNode< T_Distribution, T_Regressor, T_Splitter >
                            ( nodeId*2+1, csDim, *regressor1, Distributions[1] );
            }

            Distributions[0]=Distributions[1]=0; // the children are responsible with dealocating distribs when they don't need them anymore

            delete regressor0;
            delete regressor1;

            return moresplits;

        default:
            return false;
        }

makeleaf:
        cerr << "Something went wrong. Making node " << nodeId << "  a leaf." << endl;
        Splitter.DeleteTemporaryStatistics();

        Children[0]=Children[1]=0;

        if (buffer!=0)
        {
            delete buffer;
        }

        delete Distributions[0];
        delete Distributions[1];

        if (regressor0)
            delete regressor0;
        if (regressor1)
            delete regressor1;

        State=stable;
        return false;
    }

    /// Does the inference
    double Infer(const int* Dvars, const double* Cvars, int maxNodeId, double threshold)
    {
        if (Children[0]==0 || nodeId>maxNodeId)
        {
            // leaf node or level cut
            return Regressor.Y(Cvars+csDim);
            //return nodeId;
        }
        else
        {
            double pChild1=Splitter.ProbabilityLeft(Dvars,Cvars);

            return (pChild1>=threshold ? Children[0]->Infer(Dvars,Cvars,maxNodeId) : 0.0 )*pChild1+
                   ( 1.0-pChild1>=threshold ? Children[1]->Infer(Dvars,Cvars,maxNodeId) : 0.0 )*(1.0-pChild1);
        }
    }

    void InitializePruningStatistics(void)
    {
        pruningCost=0.0;
        pruningSamples=0;
        if (Children[0]!=0 && Children[1]!=0)
        {
            Children[0]->InitializePruningStatistics();
            Children[1]->InitializePruningStatistics();
        }
    }

    void UpdatePruningStatistics(const int* Dvars, const double* Cvars, double y /* true output */,
                                 double probability)
    {
        pruningSamples++;
        double predY=Regressor.Y(Cvars+csDim);
        pruningCost+=pow2(y-predY)*probability;
        // update pruning statistics for the proper children
        if (Children[0]!=0 && Children[1]!=0)
        {
            double pChild1=Splitter.ProbabilityFirstBranch(Dvars,Cvars);
            Children[0]->UpdatePruningStatistics(Dvars,Cvars,y,probability*pChild1);
            Children[1]->UpdatePruningStatistics(Dvars,Cvars,y,probability*(1.0-pChild1));
        }
    }

    void FinalizePruningStatistics (void)
    {
        // nothing to do
    }

    /** Returns the optimal cost for this subtree and cuts the subtree to optimal size */
    double PruneSubtree(void)
    { // double alpha /* alpha is the cost for a leaf */){
        if (Children[0] == 0 && Children[1] == 0)
        {
            // node is a leaf
            cout << "nodeID=" << nodeId << " pruningCost=" << pruningCost << endl;
            return pruningCost;
        }
        else
        {
            // node is an intermediary node
            double pruningCostChildren=Children[0]->PruneSubtree()+Children[1]->PruneSubtree();

            cout << "nodeID=" << nodeId << " pruningCost=" << pruningCost << " pruningCostChildren=" <<
            pruningCostChildren << endl;

            if (pruningCost<=pruningCostChildren)
            {
                // prune the tree here
                delete Children[0];
                Children[0]=0;
                delete Children[1];
                Children[1]=0;

                return pruningCost;
            }
            else
            {
                // tree is good as it is
                return pruningCostChildren;
            }
        }
    }

    void SaveToStream(ostream& out)
    {
        out << "{ " << nodeId << " [ ";
        if (Children[0]!=0 && Children[1]!=0)
        {
            // intermediate node
            Splitter.SaveToStream(out);
        } // else leaf, leave empty
        out << " ] ( ";
        Regressor.SaveToStream(out);
        out << " ) }";

        out << endl;

        /*
        if (parentDistribution!=0)
        parentDistribution->SaveToStream(out);
        */

        if (Children[0]!=0 && Children[1]!=0)
        {
            Children[0]->SaveToStream(out);
            Children[1]->SaveToStream(out);
        }
    }
};
}


#endif //  _CLUS_PROBABILISTICREGRESSIONTREENODE_H_

---