hipclus.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.

*/

#if !defined(_HIPCLUS_H)
#define    _HIPCLUS_H

#include    "general.h"
#include    "vec.h"
#include    "extravec.h"
#include    "fmat.h"
#include    "sphclus.h"
#include     "misc.h"
#include    "math.h"
#include    "exceptions.h"
#include    <iostream>
#include    <string>

/*    The class that implements the HiperPlan Clusters
 *    Inherit:    SphericCluster
 *    Members:
 *        C:    
 *        S:    
 *        ComputeC2s calculate the value of C2s; used when c1.. changes
 *        EigVect: the matrix of eigen vectors
 *        SmallerEigVal: the index of the smaller eigen value
 *        Elipsoidality: the ponder of the linear characteristic in the
 *            final result (takes values between 0 and 1);
 *            0 means spheric, 1 means plan
 */

namespace CLUS
{
using namespace TNT;

/// This class implements hiperclusters with only one possible output
class HiperPlanCluster:public SphericCluster
{    
protected:

    /// coefficients of the hiperplan in the form y=c0+x1*c1*..
    Vector<double> C,u;

    /// the spreading matrix; the eigen values will be found on his diagonal
    Fortran_Matrix<double> S;
    double Elipsoidality, error, slimness;
    
    void ComputeCCoef(void)
    {
        // called when u or L is modified and want to reflect modif in C
        double uy=u[inDim];
        C[0]=L[inDim];
        
        for(int j=0;j<inDim;j++)
        {
            C[j+1]=-u[j]/uy;
            C[0]+=L[j]*u[j]/uy;
        }
    }
    
    void ComputeU(void)
    {
        u[inDim]=1.0;
        L[inDim]=C[0];
        for(int j=0; j<inDim; j++)
        {
            u[j]=-C[j+1];
            L[inDim]+=C[j+1]*L[j];
        }
        Normalize(u);
    }
    
public:
    HiperPlanCluster():SphericCluster(),C(),u(),S()
    {
        Elipsoidality=0.8;
        error=0.0;
    }
    
    HiperPlanCluster(int InDim, int OutDim):SphericCluster(InDim,OutDim),C(InDim+OutDim),
            u(InDim+OutDim),S(InDim+OutDim,InDim+OutDim)
    {
        if (OutDim!=1)
            throw ErrMsg("Instances of HiperPlanClusters cannot have more than one output");
            
        Elipsoidality=0.8;
        error=0.0;
    }
    
    HiperPlanCluster(int InDim, int OutDim, int):SphericCluster(InDim,OutDim,0),C(0),
            u(InDim+OutDim),S(0,0)
    {
        if (OutDim!=1)
            throw ErrMsg("Instances of HiperPlanClusters cannot have more than one output");
        Elipsoidality=0.8;
        error=0.0;
    }
    
    void CopyFrom(HiperPlanCluster& aux)
    {
        weight=aux.weight;
        dimension=aux.dimension;
        L=aux.L;
        u=aux.u;
        error=aux.error;
        state=aux.state;
        Elipsoidality=aux.Elipsoidality;
    }

    void RandomCluster(void)
    {
        SphericCluster::RandomCluster();
        ZeroUp(u); // there is no point is choosing something else
        u[inDim]=1.0;
        ComputeCCoef();
    }
    
    void InitializeOnPoint(const double* point)
    {
        SphericCluster::InitializeOnPoint(point);
        ZeroUp(u); // there is no point is choosing something else
        u[inDim]=1.0;
        ComputeCCoef();
    }
    
    inline double InferDistance(const double* DataCache)
    {
        // the corected distance is [c1(Xj[i]-X[i])+..]^2+(Xj[i]-X[i])^2+..
        double daux=0;
        dataCache=DataCache;
        distanceInfer=0.0;
        
        for(int i=0;i<inDim;i++)
        {
            register double dif=dataCache[i]-L[i];
            distanceInfer+=pow2(dif);
            daux+=C[i+1]*dif; // corection for the inclination of the plan
        }
        
        distanceInfer+=pow2(daux);
        // we need a correction since on average the square distance to this plane
        // is error^2
        distanceInfer=(1-Elipsoidality)*distanceInfer+Elipsoidality*error;
        distanceInfer*=weight;
        
        if (distanceInfer==0.0)
            distanceInfer=10e-300;
            
        return distanceInfer;/* added with 10e-16 to avoid 0.0 distances */
    }
    
    inline double ClusDistance(const double* DataCache)
    {
        distanceClus=SphericCluster::ClusDistance(DataCache)*(1-Elipsoidality);
        double dplan=0.0;
        
        for(int j=0;j<inDim+1;j++)
            dplan+=u[j]*(dataCache[j]-L[j]);
            
        distanceClus+=Elipsoidality*pow2(dplan);//*weight;

        // sphcluster clusdistance gives the protection for 0.0 result
        return distanceClus;
    }

    /**  In the hiperplan case not only the new gravity center must be
         computed but also the dispersion matrix S
     */
    double CorrectAppartGrade(double Coef)
    {
        // we could have used spheric cluster implementation for gravity center
        // but in this way is probably slightly faster for small dimensions
        double Ai;
        if (!finite(Coef))
            Ai=1.0;
        else
            Ai=1/(distanceClus*Coef);

        double Ai2=pow2(Ai);

        assert(finite(Ai) || finite(Ai2));

        s_Ai+=Ai2;
        for (int i=0;i<inDim+1;i++)
        {
            double aux=Ai2*dataCache[i];
            s_Ai_Xj[i]+=aux;
            // fill only the upper part of S
            for(int j=i; j<inDim+1; j++)
                S(i+1,j+1)+=aux*dataCache[j];
        }
        
        return distanceClus*Ai2;
    }
    
    inline double HCorrectAppartGrade(double Coef)
    {
        double ret=CorrectAppartGrade(Coef);
        dimension+=ret;
        return ret;
    }
    
    double AdjustPrototypes(void)
    {
        double oldX, distP=0.0;
        for(int i=0;i<inDim+1;i++)
        {
            oldX=L[i];

            double aux=s_Ai_Xj[i]/s_Ai;
            L[i]=aux;
            for(int j=i; j<inDim+1; j++)
                S(i+1,j+1)-=aux*s_Ai_Xj[j]; // the corection to S
            s_Ai_Xj[i]=0.0;
            distP+=pow2(oldX-L[i]);
        }
        
        // double lambda=Min(Upper_symmetric_eigenvalue_solve(S));
        double lambda,lambda2;

        if (MinEigenvalueUpperSymmetric(S,lambda,lambda2,slimness))
        {
            if (EigenVector(S,lambda,u))
                error=lambda/s_Ai;
            else
                cerr << "OOPS could not find the eigenvector for " << lambda << endl;
            slimness/=lambda;
        }
        else
            cerr << "OOPS-eigenvalue problem not solved in " << TypeName() << endl;

        ZeroUp(S);
        s_Ai=0.0;
        ComputeCCoef(); // must stay here

        return sqrt(distP)/(inDim+1); // how much the centres moved
    }
    
    inline double HAdjustPrototypes(void)
    {
        dimension/=s_Ai;
        return AdjustPrototypes();
    }
    
    void AdjustYwithYoverD(Vector<double>& Y)
    {
        double y=C[0];
        
        for(int j=0; j<inDim; j++)
            y+=C[j+1]*dataCache[j];
            
        Y[0]+=y/distanceInfer;
    }
    
    void AdjustYwithYoverD(Vector<double>& Y, double dist)
    {
        double y=C[0];
        
        for(int j=0; j<inDim; j++)
            y+=C[j+1]*dataCache[j];
            
        Y[0]+=y/dist;
    }
    
    double ComputeLastY(void)
    {
        double y=C[0];
        
        for(int j=0; j<inDim; j++)
            y+=C[j+1]*dataCache[j];
            
        return y;
    }
    
    static string TypeName(void)
    {
        return string("HiperPlanCluster");
    }
    
    void SaveToStream(ostream& out)
    {
        ComputeCCoef();

        out << "[ ";
        // the centre of the pattern
        for(int i=0; i<inDim; i++)
            out << L[i] << " ";
        out << " ] ( ";
        // the C coefs
        for(int i=0; i<inDim+1; i++)
            out << C[i] << " ";
        out << " ) { " << Elipsoidality << " " << error
        << " " << weight << " } " << endl;
    }
    
    void LoadFromStream(istream& in)
    {
        string s;

        in >> s;
        if (s!="[")
            throw ErrMsg("[ missing in hypercluster");
        for(int i=0; i<inDim; i++)
            in >> L[i];
        in >> s;
        if (s!="]")
            throw ErrMsg("] missing in hypercluster");
        in >> s;
        if (s!="(")
            throw ErrMsg("( missing in hypercluster");
        for(int i=0; i<inDim+1; i++)
            in >> C[i];
        in >> s;
        if (s!=")")
            throw ErrMsg(") missing hypercluster");
        in >> s;
        if (s!="{")
            throw ErrMsg("{ missing in hypercluster");
        in >> Elipsoidality;
        in >> error;
        in >> weight;
        in >> s;
        if (s!="}")
            throw ErrMsg("} missing in hypercluster");
        ComputeU();
    }
    
    void SetElipsoidality(double val)
    {
        Elipsoidality=val;
    }
    
    bool SetOptionDbl(char* optName,double value)
    {
        if ( strcmp(optName,"elipsoidality")==0 )
        {
            Elipsoidality=value;
            return true;
        }
        
        return SphericCluster::SetOptionDbl(optName,value);
    }
    
    bool Split(HiperPlanCluster& right, HiperPlanCluster& save, double minDim)
    {
        if (state==Final)
            return false;
        if (state==Dead)
            throw ErrMsg("Attempt to split a dead cluster");
        if (state==Splited)
            throw ErrMsg("Attempt to split a splited cluster");

        cerr << "Slimness=" << slimness << endl;

        if (slimness <=minDim )
            return false;

        // save *this in save
        save.weight=weight;
        save.state=Splited;
        save.L=L;
        // right.C=save.C=C;
        right.u=save.u=u;
        save.error=error;
        save.dimension=dimension;
        right.Elipsoidality=save.Elipsoidality=Elipsoidality;

        state=right.state=Splitable;
        weight=right.weight=weight/2; // splited clusters are two times smaller

        Vector<double> delta(L.dim());
        GenerateRandomVector2(delta);
        // must have u*delta=0 for L-delta and L+delta to be on the plane
        double y=0.0;
        for(int j=0; j<inDim; j++)
            y+=u[j]*delta[j];
        delta[inDim]=-y/u[inDim];

        SetNormTo(delta,dimension);

        // this cluster is the left cluster
        // first the right cluster
        right.L=L+delta;
        Dif(L,L.dim(),delta);

        ComputeCCoef();
        right.ComputeCCoef();

        dimension=right.dimension=0.0;

        return true;
    }
    
    bool CheckSplit(HiperPlanCluster& left, HiperPlanCluster& right, double splitCrit)
    {
        if (error*splitCrit > left.error + right.error)
        {
            // we have a good split
            return true;
        }
        else
        {
            // we have a bad split
            right.state=Dead;
            // copy back original version
            left.weight=weight;
            left.state=Final; // we dont split this anymore
            left.L=L;
            left.u=u;
            left.error=error;
            left.Elipsoidality=Elipsoidality;
            left.dimension=0.0;
            left.ComputeCCoef();

            return false;
        }
    }
};
}


#endif    /* _HIPCLUS_H */

---