statisticsgatherers.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_STATISTICSGATHERERS_H
#define _CLUS_STATISTICSGATHERERS_H

#include "general.h"
#include "splitpointcomputation.h" // for gini and split point computation
#include "discretepermutationtransformation.h"
#include "continuouslineartransformation.h"

#include <vector>
#include <math.h>

#ifdef DEBUG_PRINT
#include <iostream>
#endif

using namespace TNT;
using namespace std;

namespace CLUS
{

enum ShiftType { labeled, unlabeled };

class IndexedValue
{
    int index;
    double value;
    
public:
    IndexedValue(int i=0, double v=0.0):
            index(i), value(v)
    {}
    
    int getIndex(void)
    {
        return index;
    }
    
    double getValue(void) const
    {
        return value;
    }
    
    bool operator < (const IndexedValue& B) const
    {
        return (value < B.getValue());
    }
    
    bool operator != (const IndexedValue& B) const
    {
        return (value != B.getValue());
    }
};

class BasicBinomialStatistics
{
protected:
    /// Counts for each value of the attribute for class 0    
    Vector<double> countsC0;
    
    /// Counts for each value of the attribute (total)
    Vector<double> counts;

    bool gainComputed;

    double ComputeSplitPoint(double N, double alpha_1)
    {
        return 0.0;
    }

public:
    BasicBinomialStatistics(int DomainSize=0):
            countsC0(DomainSize), counts(DomainSize), gainComputed(false)
    {}

    double getCount()
    {
        double total=0;
        for (int i=0; i<countsC0.size(); i++)
            total+= counts[i];
        return total;
    }

    void Print(void)
    {
        cout << "CountC0 [ ";
        for (int i=0; i<countsC0.size(); i++)
            cout << countsC0[i] << " ";
        cout << " ]" << endl;

        cout << "Count [";
        for (int i=0; i<counts.size(); i++)
            cout << counts[i] << " ";
        cout << " ]" << endl;
    }

    void ResetDomainSize(int DomainSize)
    {
        countsC0.newsize(DomainSize);
        counts.newsize(DomainSize);
        for (int i =0; i< DomainSize; i++)
        {
            countsC0[i]=0.0;
            counts[i]=0.0;
        }
    }

    void UpdateStatistics(int value, int classLabel, double probability=1.0)
    {
        assert(probability>=0.0 && probability<=1.0);

        if (classLabel==0)
            countsC0[value]+=probability;
        counts[value]+=probability;
#ifdef DEBUG_PRINT

        cout << " value=" << value <<  " counts[value]=" << counts[value]
        << " countsC0[value]=" << countsC0[value] << endl;
#endif

        assert (isfinite(counts[value]));
    }

    void UpdateStatisticsP(int value, double p0, double p1)
    {
        countsC0[value]+=p0;
        counts[value]+=p0+p1;
        assert (isfinite(counts[value]));
    }
};


class BinomialStatistics : public BasicBinomialStatistics
{
protected:
    /// the values in the left set; split point
    Vector<int> Split; 

    double ComputeSplitPoint(double N, double alpha_1)
    {
        return DiscreteGiniGain(countsC0, counts, N, alpha_1, Split);
    }

public:
    BinomialStatistics(int DomainSize=0): BasicBinomialStatistics(DomainSize),
            Split(0)
    {}

    Vector<int>& GetSplit(void)
    {
        if (gainComputed)
            return Split;
        else
        {
            ComputeGiniGain();
            return Split;
        }
    }

    double ComputeGiniGain(void)
    {
        // create a temporary vector of doubles out of countsC0 so that
        // DiscreteGiniGain function can be used
        gainComputed = true;
        double N=0.0; // keeps the total count
        double NC0=0.0;
        for (int i=0; i<counts.dim(); i++)
        {
            N+=counts[i];
            NC0+=countsC0[i];
        }
#ifdef DEBUG_PRINT
        cout << "Discrete. ComputeGiniGain " << N << " " << NC0 << "\t";
#endif

        double alpha_1=NC0/N;
        //cerr << alpha_1 << "?" << endl;
        if (alpha_1!=alpha_1 || N==0.0 )
        {
            //cerr << "no discrete data"<< endl;
            return 0;
        }

        double val = DiscreteGiniGain(countsC0, counts, N, alpha_1, Split);
#ifdef DEBUG_PRINT

        cout << val << endl;
#endif

        return val;
    }
    
    // entries r_i/q_i
    vector<IndexedValue> getConditionalProbs(void)
    {
        vector<IndexedValue> myVals(counts.size());
        double C0_total=0;//how many C0 examples were seen
        double total = 0;//how many examples were seen
        for (int i=0; i<counts.size(); i++)
        {
            C0_total+= countsC0[i];
            total += counts[i];
        }
        // cerr << "0's: " << C0_total << " total: " << total << endl;
        for (int i=0; i<counts.size(); i++)
        {
            double val = (double)(countsC0[i]/C0_total)/(counts[i]/total);
            //  cerr << "(" << i << ", " << val << endl;
            IndexedValue iv(i, val);
            myVals[i] = iv;
        }
        return myVals;
    }

    Permutation ComputeShift(ShiftType shift, BinomialStatistics& aux)
    {
        // compute the permutation such that the rank of the
        // probabilities to see classLabel=0 in both is the same
        if(!gainComputed)
        {
            ComputeGiniGain();
        }
        vector<IndexedValue> myVals = getConditionalProbs();
        vector<IndexedValue> refVals = aux.getConditionalProbs();
        sort(myVals.begin(), myVals.end());
        sort(refVals.begin(), refVals.end());
        Vector<int> perm(counts.size());
        for (int i =0; i<counts.size(); i++)
        {
            int myNextIndex = myVals[i].getIndex();
            int refNextIndex = refVals[i].getIndex();
            //  cerr << "mapping " << refNextIndex << " to " <<  myNextIndex << endl;
            perm[refNextIndex] = myNextIndex;
        }
        return Permutation(perm);
    }

    void AddWeightedStatistics(BinomialStatistics& aux, double weight)
    {
        assert(counts.dim() == aux.counts.dim());

        if (weight==0.0)
            return;

        for (int i=0; i<counts.dim(); i++)
        {
            counts[i]+=weight*aux.counts[i];
            countsC0[i]+=weight*aux.countsC0[i];
        }
    }

    void CorrectWeightedStatistics(double totalWeight)
    {
        for (int i=0; i<counts.dim(); i++)
        {
            counts[i]/=totalWeight;
            countsC0[i]/=totalWeight;
        }
    }

    void AddStatisticsShifted(BinomialStatistics& aux, Permutation shift)
    {
        assert(counts.dim() == aux.counts.dim());
        // apply the shift when adding the statistics
        for (int i=0; i<counts.dim(); i++)
        {
            int shiftedIndex=shift.ApplyPermutation(i);
            counts[i]+=aux.counts[shiftedIndex];
            countsC0[i]+=aux.countsC0[shiftedIndex];
        }
    }

};

class ProbabilisticBinomialStatistics : public BasicBinomialStatistics
{
protected:
    Vector<double> probSet; // probability to be at the left for each value

public:
    ProbabilisticBinomialStatistics(int DomainSize=0): BasicBinomialStatistics(DomainSize),
            probSet(DomainSize)
    {}

    Vector<double>& GetProbabilitySet(void)
    {
        if (gainComputed)
            return probSet;
        else
        {
            ComputeGiniGain();
            return probSet;
        }
    }

    double ComputeGiniGain(void)
    {
        // create a temporary vector of doubles out of countsC0 so that
        // DiscreteGiniGain function can be used
        gainComputed = true;
        double N=0.0; // keeps the total count
        double NC0=0.0;
        for (int i=0; i<counts.dim(); i++)
        {
#ifdef DEBUG_PRINT
            cout << " counts[" << i << "]=" << counts[i]
            << " countsC0[" << i << "]=" << countsC0[i] << endl;
#endif

            N+=counts[i];
            NC0+=countsC0[i];
        }
#ifdef DEBUG_PRINT
        cout << "Discrete. ComputeGiniGain " << N << " " << NC0 << "\t";
#endif

        double alpha_1=NC0/N;
        //cerr << alpha_1 << "?" << endl;
        if (alpha_1!=alpha_1 || N==0.0 )
        {
            //cerr << "no discrete data"<< endl;
            return 0;
        }

        double val = ProbabilisticDiscreteGiniGain(countsC0, counts, N, alpha_1, probSet);
#ifdef DEBUG_PRINT

        cout << val << endl;
#endif

        return val;
    }


};

class NormalStatistics
{
protected:
    double countC0;
    double countC1;

    /// maintain the sum of values for class label 0
    double sumC0; 
    double sumC1;

    /// maintain the sum of squares for class label 0
    double sum2C0; 
    double sum2C1;
    double split;

    /// Solution returned by QDA
    int whichSol; 
    double splitVariance;

    bool splitVarComputed;



public:
    NormalStatistics(int dummy=0)
    {
        Reset();
    }

    void Reset(void)
    {
        countC0=0.0;
        countC1=0.0;
        sumC0=0.0;
        sumC1=0.0;
        sum2C0=0.0;
        sum2C1=0.0;
        splitVarComputed=false;
    }

    void Print(void)
    {
        cout << "Statistics: " << countC0 << " " << countC1 << " ";
        cout << sumC0 << " " << sumC1 << " " << sum2C0 << " " << sum2C1 << endl;
    }

    double getVariance(void)
    {
        return (sum2C0+sum2C1)/(countC0+countC1) - pow2((sumC0+sumC1)/(countC0+countC1));
    }

    double getSplitVariance(void)
    {
        if (!splitVarComputed)
        {
            ComputeGiniGain();
        }
        //return 1.0e-32; // FIXTHIS
        return splitVariance;// is this variance???
    }

    void UpdateStatistics(double value, int classLabel, double probability=1.0)
    {
        if (classLabel==0)
        {
            countC0+=probability;
            sumC0+=probability*value;
            sum2C0+=probability*pow2(value);
        }
        else
        {
            countC1+=probability;
            sumC1+=probability*value;
            sum2C1+=probability*pow2(value);
        }
        assert(isfinite(sumC0));// NaN
    }

    double GetCenter(void)
    {
        return (sumC0+sumC1)/(countC0+countC1);
    }

    double ComputeGiniGain(void)
    {
        if (countC0+countC1==0)
            return 0.0;

        double alpha_1=(double)countC0/(countC0+countC1);
        double alpha_2=(double)countC1/(countC0+countC1);
        double eta1=sumC0/countC0;
        double eta2=sumC1/countC1;
        double var1=sum2C0/countC0-pow2(eta1);
        double var2=sum2C1/countC1-pow2(eta2);

#ifdef DEBUG_PRINT

        cout << "ComputeGiniGain " << alpha_1 << " " << alpha_2 << " " << eta1 << " ";
        cout << eta2 << " " << var1 << " " << var2 << "\t";
#endif

        split=UnidimensionalQDA(alpha_1,eta1,var1,alpha_2,eta2,var2, whichSol);
        assert(isfinite(split));
        splitVariance=UnidimensionalQDAVariance(countC0, eta1, var1,
                                                countC1, eta2, var2, whichSol);
        splitVarComputed=true;

        // compute gini like in function ComputeSeparatingHyperplane_Anova
        double p11=alpha_1*PValueNormalDistribution(eta1, sqrt(var1), split);
        double p1_=p11+alpha_2*PValueNormalDistribution(eta2, sqrt(var2), split);
        double gini=BinaryGiniGain(p11,alpha_1,p1_);
#ifdef DEBUG_PRINT

        cout << split << " " << gini << endl;
#endif

        return gini;
    }

    double GetSplit(void)
    {
        return split;
    }

    /* --------------------- Stuff only for MultiDecTrees -------------*/

    bool nonZero(void)
    {
        return (countC0!=0 && countC1!=0);
    }

    bool hasData(void)
    {
        return (countC0!=0 && countC1!=0);
    }

    double ComputeShift(ShiftType shift, NormalStatistics& aux)
    {
        //      cerr << "computing shift" << endl;
        if (shift==unlabeled)
            return ( (aux.sumC0+aux.sumC1)/(aux.countC0+aux.countC1) -
                     (sumC0+sumC1)/(countC0+countC1) );
        /* The ComputeGiniGain for (*this) has to be called before */
        if (shift==labeled)
        {
            aux.ComputeGiniGain();
            assert (!(split!=split));
            return aux.split-split;
        }
        return 0.0;
    }

    void AddStatisticsShifted(NormalStatistics& aux, double shift)
    {
        // the shift has to be substracted from the statistics of aux
        countC0+=aux.countC0;
        countC1+=aux.countC1;
        // simulate the situation where the shift is substracted from every datapoint
        sumC0+=aux.sumC0-aux.countC0*shift;
        assert(!(sumC0!=sumC0)); //Nan
        assert(!(shift!=shift));
        sumC1+=aux.sumC1-aux.countC1*shift;
        sum2C0+=aux.sum2C0-2*aux.sumC0*shift+pow2(shift)*aux.countC0;
        sum2C1+=aux.sum2C1-2*aux.sumC1*shift+pow2(shift)*aux.countC1;
    }

    void AddWeightedStatistics(NormalStatistics& aux, double weight)
    {
        if (weight==0.0)
            return;

        countC0+=aux.countC0*weight;
        countC1+=aux.countC1*weight;
        sumC0+=aux.sumC0*weight;
        sumC1+=aux.sumC1*weight;
        sum2C0+=pow2(weight)*(aux.sum2C0*aux.countC0-pow2(aux.sumC0));
        sum2C1+=pow2(weight)*(aux.sum2C1*aux.countC1-pow2(aux.sumC1));
    }

    void CorrectWeightedStatistics(double totalWeight)
    {
        countC0/=totalWeight;
        countC1/=totalWeight;
        sumC0*=countC0/totalWeight;
        sumC1*=countC1/totalWeight;
        sum2C0=sum2C0/pow2(countC0*totalWeight)+pow2(sumC0/countC0);
        sum2C1=sum2C1/pow2(countC1*totalWeight)+pow2(sumC1/countC1);
    }

    bool negMeanLessThanPos(void)
    {
        double negMean = sumC0/countC0;
        double posMean = sumC1/countC1;
        return (negMean < posMean);
    }

    bool labeledMeansSignificant(void)
    {
        double negMean = sumC0/countC0;
        double posMean = sumC1/countC1;
        double diff = posMean-negMean;
        if (!negMeanLessThanPos())
        {
            diff = -1.0*diff;
        }

        /*
        double posVar = getLCVariance(true);
        double negVar = getLCVariance(false);
        */
        if (diff >= 0)// .0000000001*(posVar+negVar))
        {
            return true;
        }
        else
            return false;

    }

    double getLabeledCenter(bool b)
    {
        if (b)
            return sumC1/countC1;
        else
            return sumC0/countC0;
    }

    double  getLCVariance(bool b)
    {
        if (b)
        {
            return sum2C1/countC1 - pow2(sumC1/countC1);
        }
        else
            return sum2C0/countC0 - pow2(sumC0/countC0);
    }

    double getcountC0()
    {
        return countC0;
    }
    
    double getcountC1()
    {
        return countC1;
    }

};


/** Class implements a multidimentional normal distribution
 */
class MultidimNormalStatistics
{
protected:
    int dim;
    double s_p1, s_p2; // Sum p1i, Sum p2i
    Vector<double> s_p1_x, s_p2_x; // Sum p1i*x, Sum p2i*x
    Fortran_Matrix<double> S_p1_xxT, S_p2_xxT; // Sum p1i*xx^T, Sum p2i*xx^T

    int SplitVariable;
    Vector<double> SeparatingHyperplane;

public:
    MultidimNormalStatistics()
    {
        dim=0;
    } // all the default constructors called for vectors and matrices

    enum SeparationType { ANOVA=0, LDA=1, QDA=2 };

    double GetS_P1(void)
    {
        return s_p1;
    }
    
    double GetS_P2(void)
    {
        return s_p2;
    }

    /** changes the dimention and reinitializes everything */
    void Resize(int Dim)
    {
        if (dim!=Dim)
        {
            dim=Dim;

            s_p1_x.newsize(dim);
            s_p2_x.newsize(dim);
            S_p1_xxT.newsize(dim,dim);
            S_p2_xxT.newsize(dim,dim);
            SeparatingHyperplane.newsize(dim+1);
        }

        if (dim>0) /* otherwise somebody is just freeing the memory */
            Reset();
    }

    void Reset(void)
    {
        s_p1=s_p2=0.0;
        s_p1_x=0.0;
        s_p2_x=0.0;
        S_p1_xxT=0.0;
        S_p2_xxT=0.0;
    }

    void UpdateStatistics(const double* values, double p1, double p2)
    {
        if (!finite(p1+p2))
            return;

        s_p1+=p1;
        s_p2+=p2;

        for(int i=0; i<dim; i++)
        {
            double aux1=p1*values[i];
            double aux2=p2*values[i];
            s_p1_x[i]+=aux1;
            s_p2_x[i]+=aux2;
            for(int j=0; j<dim; j++)
            {
                S_p1_xxT(i+1,j+1)+=aux1*values[j];
                S_p2_xxT(i+1,j+1)+=aux2*values[j];
            }
        }
    }

    /* must be called before the split point can be computed */
    void UpdateParameters(void)
    {
        for(int i=0; i<dim; i++)
        {
            // compute mu1, mu2 in s_p1_x and s_p2_x
            s_p1_x[i]/=s_p1;
            s_p2_x[i]/=s_p2;
        }
        // compute Sigma1 and Sigma2 in S_p1_xxT and S_p2_xxT
        for(int i=0; i<dim; i++)
        {
            for(int j=0; j<dim; j++)
            { // AlinTODO: scan only half of the matrix
                S_p1_xxT(i+1,j+1)=S_p1_xxT(i+1,j+1)/s_p1-s_p1_x[i]*s_p1_x[j];
                S_p2_xxT(i+1,j+1)=S_p2_xxT(i+1,j+1)/s_p2-s_p2_x[i]*s_p2_x[j];

                // fix in S_xxT the almost 0.0 entries. This helpes a lot
                if (fabs(S_p1_xxT(i+1,j+1))<SMALL_POZ_VALUE)
                {
                    S_p1_xxT(i+1,j+1)=0.0;
                }
                if (fabs(S_p2_xxT(i+1,j+1))<SMALL_POZ_VALUE)
                {
                    S_p2_xxT(i+1,j+1)=0.0;
                }

            }
        }
    }

    int GetSplitVariable(void)
    {
        return SplitVariable;
    }

    Vector<double>& GetSeparatingHyperplane(void)
    {
        return SeparatingHyperplane;
    }

    double ComputeGiniGain(int type = LDA)
    {
        double gini, alpha_1, alpha_2;

        if (dim<=1)
            return 0.0;

        SplitVariable=-1;

        if (s_p1+s_p2==0.00)
            return 0.0;

        // compute the best oblique split
        alpha_1 =s_p1/(s_p1+s_p2);
        alpha_2 =s_p2/(s_p1+s_p2);
        double mass=s_p1+s_p2;

        UpdateParameters();

        /* All the SeparatingHyperplane computation funcions are normalized
        with respect to variance of the split point */

        switch (type)
        {
        case ANOVA:
            gini = ComputeSeparatingHyperplane_Anova(mass,
                    alpha_1, s_p1_x, S_p1_xxT,
                    alpha_2, s_p2_x, S_p2_xxT, SeparatingHyperplane);
            break;
        case LDA:
            gini = ComputeSeparatingHyperplane_LDA(mass,
                                                   alpha_1, s_p1_x, S_p1_xxT,
                                                   alpha_2, s_p2_x, S_p2_xxT, SeparatingHyperplane);
            break;
        case QDA:
            gini = ComputeSeparatingHyperplane_QDA(mass,
                                                   alpha_1, s_p1_x, S_p1_xxT,
                                                   alpha_2, s_p2_x, S_p2_xxT, SeparatingHyperplane);
            break;
        default:
            gini = 0.0; // should never be reached
        }

        SplitVariable=-1;
        {
            // See if we have a simple split
            bool isSimple = true;
            int posSplitVar = -1;
            for (int i=0; i<SeparatingHyperplane.dim()-1; i++)
            {
                double cval=SeparatingHyperplane[i+1];
                if (cval != 0.0)
                    if (posSplitVar==-1)
                        posSplitVar=i;
                    else
                    {
                        isSimple=false;
                        break;
                    }
                else
                    continue;
            }

            if (isSimple && posSplitVar!=-1)
            {
                SplitVariable = -(posSplitVar+2);
                cout << "Split variable is simple=" << posSplitVar << endl;
            }
            else
            {
                cout << "Simple separation detection failed at: " << posSplitVar << endl;
            }

        }

        return gini;
    }

    double MaxGini(void)
    {
        return 2*s_p1*s_p2/pow2(s_p1+s_p2);
    }

};

}

#endif // _CLUS_STATISTICSGATHERERS_H

---