cfentry.c

数据挖掘经典的hierarchial clustering algorithm

/****************************************************************
File Name:   cfentry.C
Author: Tian Zhang, CS Dept., Univ. of Wisconsin-Madison, 1995

               Copyright(c) 1995 by Tian Zhang

                   All Rights Reserved

Permission to use, copy and modify this software must be granted
by the author and provided that the above copyright notice appear
in all relevant copies and that both that copyright notice and this
permission notice appear in all relevant supporting documentations.

Comments and additions may be sent the author at [e-mail].

******************************************************************/

#include "global.h"
#include "util.h"
#include "vector.h"
#include "rectangle.h"
#include "cfentry.h"

Entry::Entry() {}

void Entry::Init(short d) {

n=0; sx.Init(d); sxx=0;

#ifdef RECTANGLE
rect.Init(d);
#endif RECTANGLE
}

void Entry::Reset() {

n=0; sx.Reset(); sxx=0;

#ifdef RECTANGLE
rect.Reset();
#endif RECTANGLE
}

Entry::~Entry() {}

// correct copy constructor
Entry::Entry(const Entry& r) {
if (this!=&r) {
        n = r.n;
        sx = r.sx;
        sxx = r.sxx;

#ifdef RECTANGLE
        rect = r.rect;
#endif RECTANGLE
}
}

// correct assignment operator
void Entry::operator=(const Entry& r) {
if (this!=&r) {
        n = r.n;
        sx = r.sx;
        sxx = r.sxx;

#ifdef RECTANGLE
        rect = r.rect;
#endif RECTANGLE
        }
}

void Entry::operator=(const int val) {
        n = val;
        sx = val;
        sxx = val;

#ifdef RECTANGLE
        rect = val;
#endif RECTANGLE
        }

void Entry::operator=(const Vector &v) {
        n = 1;
        sx = v;
        sxx = sx && sx;

#ifdef RECTANGLE
        rect = v;
#endif RECTANGLE
        }

short Entry::Dim() const {return sx.Dim();}

int Entry::N() const {return n;}

void Entry::SX(Vector& tmpsx) const {tmpsx = sx;}

double Entry::SXX() const {return sxx;}

void Entry::X0(Vector &tmpx0) const {tmpx0.Div(sx,n);}

#ifdef RECTANGLE
void Entry::Rect(Rectangle& tmprect) const {tmprect = rect;}
#endif RECTANGLE

// Assume:
// K(t) is standard normal distribution.
// gi(x) is normal distribution with its own xi and ri.
// wi(x) expanded only to o(h^2) with Taylor Expansions.
// work for $d$ dimensions.

double Entry::Norm_Kernel_Density_Effect(const Vector &x, double h) const {
double ri;
double tmp0=0;
double tmp1=0;

// w(x) = Norm(x0, h^2+r^2)

ri = sqrt(this->Radius()/sx.dim+h*h);

for (short i=0; i<sx.dim; i++) {
    tmp0 = (x.value[i]-sx.value[i]/n)/ri;
    tmp1 -= tmp0*tmp0/2.0;
    }
return exp(tmp1)/pow(sqrt(2*PI)*ri,sx.dim);
}

// work for 1 dimension, need to generalize
double Entry::Norm_Kernel_Prob_Effect(const Vector& a, double h) const
{
double r2=Radius();
return 0.5*(1.0+erf((a.value[0]-sx.value[0]/n)/(sqrt(2.0*(r2+h*h)))));
}

// Works for 1-dim presently, need to generalize to d-dim.
double Norm_Kernel_Rf_Effect(const Entry& ent1, const Entry& ent2, double h)
{
double ri,rj;
double tmp0,tmp1,tmp2;

ri = ent1.Radius();
rj = ent2.Radius();
tmp0 = 2.0*h*h+ri+rj;
tmp1 = ent1.sx.value[0]/ent1.n-ent2.sx.value[0]/ent2.n;
tmp2 = tmp1*tmp1;
tmp1 = tmp2/tmp0;
return 1.0/(sqrt(2.0*PI*pow(tmp0,5.0)))*exp(-tmp1/2.0)*((tmp1-3.0)*(tmp1-3.0)-6);
}

// Assume:
// K(t) is standard normal distribution.
// gi(x) is uniform distribution with its own xi and ri.
// wi(x) is actually integrated.
// work for $d$ dimensions.

double Entry::Unif_Kernel_Density_Effect(const Vector &x, double h) const {
double ri;
double tmp0;
double tmp1;
double tmp=1.0;

ri = this->Radius()/sx.dim;

if (ri<=0) return this->Norm_Kernel_Density_Effect(x,h);

else {  ri = sqrt(ri);
        tmp0 = sqrt(3.0)*ri;
        tmp1 = sqrt(2.0)*h;
        for (short i=0; i<sx.dim; i++) {
                tmp *=(erf((sx.value[i]/n+tmp0-x.value[i])/tmp1)
                      -erf((sx.value[i]/n-tmp0-x.value[i])/tmp1))
                      /(4.0*tmp0);
                }
        return tmp;
        }
}

// work for 1 dimension, need to generalize
double Entry::Unif_Kernel_Prob_Effect(const Vector &a, double h) const
{
double tmp0, tmp1, tmp2, tmp3, tmp4;
double r2=Radius();
if (r2<=0)
  return 0.5*(1.0+erf((a.value[0]-sx.value[0]/n)/(sqrt(2.0)*h)));
else {
  tmp0=sqrt(3.0*r2);
  tmp1=sqrt(2.0)*h;
  tmp2=a.value[0]-sx.value[0]/n;
  tmp3=tmp2-tmp0;
  tmp4=tmp2+tmp0;
  return 0.5+1.0/(4.0*tmp0)*(-tmp3*erf(tmp3/tmp1)+tmp4*erf(tmp4/tmp1)
                             +sqrt(2.0/PI)*h*(-exp(-tmp3/tmp1*tmp3/tmp1)
                                              +exp(-tmp4/tmp1*tmp4/tmp1)));
  }
}

// Works for 1-dim presently, need to generalize to d-dim.
double Unif_Kernel_Rf_Effect(const Entry& ent1, const Entry& ent2, double h)
{
double ri,rj;
double xi,xj;
double ai,bi,aj,bj,tmp0,tmp1,tmp2;

ri=sqrt(ent1.Radius());
rj=sqrt(ent2.Radius());
xi=ent1.sx.value[0]/ent1.n;
xj=ent2.sx.value[0]/ent2.n;

tmp0 = sqrt(3.0)*ri;
ai=xi-tmp0;
bi=xi+tmp0;
tmp1 = sqrt(3.0)*rj;
aj=xj-tmp1;
bj=xj+tmp1;

if (ri<=0 && rj<=0) {
        tmp0=2*h*h;
        tmp1=(xi-xj)*(xi-xj);
        return 1.0/sqrt(2.0*PI*pow(tmp0,5.0))*exp(-tmp1/(2.0*tmp0))*((tmp1/tmp0-3.0)*(tmp1/tmp0-3.0)-6.0);
        }

if (ri>0 && rj>0) {
        tmp2=0.0;
        tmp0=h*h;
        tmp1=(bi-bj)*(bi-bj);
        tmp2+=exp(-tmp1/(4.0*tmp0))*(tmp0-tmp1/2.0);
        tmp1=(ai-bj)*(ai-bj);
        tmp2-=exp(-tmp1/(4.0*tmp0))*(tmp0-tmp1/2.0);
        tmp1=(bi-aj)*(bi-aj);
        tmp2-=exp(-tmp1/(4.0*tmp0))*(tmp0-tmp1/2.0);
        tmp1=(ai-aj)*(ai-aj);
        tmp2+=exp(-tmp1/(4.0*tmp0))*(tmp0-tmp1/2.0);
        return tmp2/(48*sqrt(PI)*ri*rj*h*h*h*h*h);
        }

if (ri<=0 && rj>0) {
        tmp2=0.0;
        tmp0=h*h;
        tmp1=(xi-bj)/2.0;
        tmp2+=exp(-1.0/tmp0*tmp1*tmp1)*(1.5*tmp1-tmp1*tmp1*tmp1/tmp0);
        tmp1=(xi-aj)/2.0;
        tmp2-=exp(-1.0/tmp0*tmp1*tmp1)*(1.5*tmp1-tmp1*tmp1*tmp1/tmp0);
        return tmp2/(-4.0*sqrt(3.0*PI)*rj*h*h*h*h*h);
        }

if (ri>0 && rj<=0) {
        tmp2=0.0;
        tmp0=h*h;
        tmp1=(xj-bi)/2.0;
        tmp2+=exp(-1.0/tmp0*tmp1*tmp1)*(1.5*tmp1-tmp1*tmp1*tmp1/tmp0);
        tmp1=(xj-ai)/2.0;
        tmp2-=exp(-1.0/tmp0*tmp1*tmp1)*(1.5*tmp1-tmp1*tmp1*tmp1/tmp0);
        return tmp2/(-4.0*sqrt(3.0*PI)*ri*h*h*h*h*h);
        }
}

double Entry::Diameter() const {
        if (n<=1) return 0;
        else // return 2*(n*sxx-(sx&&sx))/(n*(n-1));
             // for performance reason
             {double tmp1,tmp0=0;
              for (short i=0; i<sx.dim; i++)
                tmp0 += sx.value[i]/n*sx.value[i]/(n-1);
              tmp1=2*(sxx/(n-1.0)-tmp0);
              if (tmp1<=0.0) return 0.0;
              else return tmp1;
              }
}

double Entry::Radius() const {
        if (n<=1) return 0;
        else // return(sxx/n-((sx/n)&&(sx/n)));
             // for performance reason
                {double tmp0, tmp1=0;
                 for (short i=0; i<sx.dim; i++) {
                        tmp0 = sx.value[i]/n;
                        tmp1 += tmp0*tmp0;
                        }
                 tmp0=sxx/n-tmp1;
                 if (tmp0<=0.0) return 0.0;
                 else return tmp0;
                }
}

double Entry::Fitness(short ftype) const {
        switch (ftype) {
        case AVG_DIAMETER: return Diameter();
        case AVG_RADIUS:   return Radius();
        default: print_error("Entry::Fitness", "Invalid fitness type");
        }
}

void Entry::operator+=(const Entry &ent) {
        n += ent.n;
        sx += ent.sx;
        sxx += ent.sxx;

#ifdef RECTANGLE
        rect += ent.rect;