|| || || || || || || ||

---

MatKMeans.h

Go to the documentation of this file.

#ifndef Impala_Core_Matrix_MatKMeans_h
#define Impala_Core_Matrix_MatKMeans_h

#include "Core/Matrix/MatFunc.h"
#include "Core/Matrix/Vecs.h"

namespace Impala
{
namespace Core
{
namespace Matrix
{


// calculate vec distance
inline double
my_vec_distance(double* v1, double* v2, int nElem) 
{
    double s=0,tmp ;
    for (int i=0;i<nElem;i++) 
        s += (tmp = v1[i] - v2[i],tmp*tmp);
    return sqrt(s);
}

// calculate vec distance, given a current minimum
inline double
my_vec_distance(double* v1, double* v2, int nElem, double currentMin) 
{
    if (currentMin>0)
    {
        double currentMinSquared = currentMin*currentMin;
        double s=0,tmp ;
        int i=0 ;
        while (s<currentMinSquared && i<nElem)
        {
                s += (tmp = v1[i] - v2[i],tmp*tmp);
                i++;
        }
        return i==nElem?sqrt(s):currentMin;
    }
    return currentMin;  
}

template<class ArrayT>
void
MatKMeans(ArrayT* m, ArrayT* c, VecScalarReal64* imap,
          int numloop, bool cIsInitialized, int& clustersFound,
          VecScalarInt32* vectorsInC)
{
    typedef typename ArrayT::StorType StorT;
    typedef typename ArrayT::ArithType ArithT;

    int nrVectors = MatNrRow(m);
    int nrDimensions = MatNrCol(m);
    int nrClusters = MatNrRow(c); 

    int i,j;
    StorT* cPtr;
    StorT* mPtr;
    double* imapPtr = VecE(imap, 0);

    if (!cIsInitialized)
    {
        // the cluster matrix is not initialized, initialize it with some vectors.
        int steps = nrVectors / nrClusters;
        cPtr = MatE(c, 0, 0);
        for (i=0 ; i<nrClusters ; i++)
        {
            mPtr = MatE(m, i*steps, 0);
            for (j=0 ; j<nrDimensions ; j++)
            {
                *cPtr++ = *mPtr++;
            }
        }
    } 

    bool contin = true;
    ArithT err = 0, errOld = 0;
    int emptyclust = 0;
    StorT* nextC = new StorT[nrClusters*nrDimensions];
    int* nrInNextC = new int[nrClusters];
    StorT* nextPtr;
    while (contin)
    {
        emptyclust = 0;
        err = 0;
        contin = false;
        nextPtr = nextC;
        for (i=0 ; i<nrClusters ; i++)
        {
            nrInNextC[i] = 0;
            for (j=0 ; j<nrDimensions ; j++)
                *nextPtr++ = 0;
        }

        for (i=0 ; i<nrVectors ; i++)
        {
            int myClus = 0;
            ArithT d = 0;
            mPtr = MatE(m, i, 0);
            ArithT dm = my_vec_distance(mPtr, MatE(c, 0, 0), nrDimensions);
            for (j=1 ; j<nrClusters ; j++)
            {
                //d = vec_distance(mPtr, MatE(c, j, 0), nrDimensions, dm);
                d = my_vec_distance(mPtr, MatE(c, j, 0), nrDimensions);
                if (d<dm)
                {
                    dm = d;
                    myClus = j;
                }
            }
            imapPtr[i] = myClus;
            err += dm;

            // store this vector's information in nextC for next iteration
            nrInNextC[myClus] += 1;
            nextPtr = nextC + myClus*nrDimensions;
            for (j=0 ; j<nrDimensions ; j++)
                *nextPtr++ += *mPtr++;
        }
        
        numloop--;
        // recalculate the position of the cluster vectors
        cPtr = MatE(c, 0, 0);
        nextPtr = nextC;
        for (i=0 ; i<nrClusters ; i++)
        {
            if (nrInNextC[i] > 0)
            {
                for (j=0 ; j<nrDimensions ; j++)
                    *cPtr++ = *nextPtr++ / nrInNextC[i];
            }
            else
            {
                for (j=0 ; j<nrDimensions ; j++)
                    *cPtr++ = *nextPtr++;
                emptyclust++ ;
                printf("empty cluster: %d\n",i);
            }
        }

        //printf("%d, ",numloop);
        /*
        cPtr = c;
        printf("\nloop %d\n", numloop);
        for (i=0 ; i<nrClusters ; i++)
        {
            printf("c %d : ", i);
            for (j=0 ; j<nrDimensions ; j++)
                printf("%f ", *cPtr++);
            printf("\n");
        }
        */
        // quit if the numLoops is reached or if the error stays constant
        if (err != errOld && numloop > 0)
            contin = true;

        errOld = err;
    } // end while(contin)
    //printf("\n");

    for (i=0 ; i<nrClusters ; i++)
        //vectorsInC[i] = nrInNextC[i];
        *VecE(vectorsInC, i) = nrInNextC[i];
    clustersFound = nrClusters - emptyclust;

    delete nextC;
    delete nrInNextC;
}

} // namespace Matrix
} // namespace Core
} // namespace Impala

#endif

---