SVClustering.cpp

Go to the documentation of this file.

#include "SVClustering.h"
#include "L2Distance.h"

namespace damina
{

        SVClustering::SVClustering(DataSet *train) {
                this->eucDist = new L2Distance();
                this->trainer = new OneClassSVM();
                this->trainer->setKernel(KernelType::gaussian);
                this->trainer->setCacheSize(100);
                this->trainer->setTolerance(0.000001);
                this->trainer->setTrainingSet(train);
                this->trainer->setKernelWidth(initialKernelWidth());
                this->quadratic = NULL;
                this->clusterizeBSV = true;
        }
        
        SVClustering::SVClustering(struct svm_problem *p) {
                this->eucDist = new L2Distance();
                this->trainer = new OneClassSVM();
                this->trainer->setKernel(KernelType::gaussian);
                this->trainer->setCacheSize(100);
                this->trainer->setTolerance(0.000001);
                this->trainer->setProblem(p);
                this->trainer->setKernelWidth(initialKernelWidth());
                this->quadratic = NULL;
                this->clusterizeBSV = true;
        }
        
        SVClustering::SVClustering(double C, DataSet *train) {
                this->eucDist = new L2Distance();
                this->trainer = new OneClassSVM();
                this->trainer->setKernel(KernelType::gaussian);
                this->trainer->setCacheSize(100);
                this->trainer->setTolerance(0.000001);
                
                this->trainer->setTrainingSet(train);
                
                this->trainer->setKernelWidth(initialKernelWidth());

                
                this->trainer->setNU((1 / (train->getSize() * C)));
                this->C = C;
                this->quadratic = NULL;
                this->clusterizeBSV = true;
        }
        
        SVClustering::SVClustering(double C, struct svm_problem *prob) {
                this->eucDist = new L2Distance();
                this->trainer = new OneClassSVM();
                this->trainer->setKernel(KernelType::gaussian);
                this->trainer->setCacheSize(100);
                this->trainer->setTolerance(0.000001);
                
                this->trainer->setProblem(prob);
                
                this->trainer->setKernelWidth(initialKernelWidth());
                
                this->trainer->setNU((1 / (prob->l * C)));
                this->C = C;
                this->quadratic = NULL;
                this->clusterizeBSV = true;
        }
        
        SVClustering::SVClustering(double C, double w, DataSet *train) {
                this->eucDist = new L2Distance();
                this->trainer = new OneClassSVM();
                this->trainer->setKernel(KernelType::gaussian);
                this->trainer->setCacheSize(100);
                this->trainer->setTolerance(0.000001);
                
                this->trainer->setTrainingSet(train);
                
                this->trainer->setNU((1 / (train->getSize() * C)));
                this->C = C;
                
                this->trainer->setKernelWidth(w);
                this->quadratic = NULL;
                this->clusterizeBSV = true;
        }
        
        SVClustering::SVClustering(double C, double w, struct svm_problem *prob) {
                this->eucDist = new L2Distance();
                this->trainer = new OneClassSVM();
                this->trainer->setKernel(KernelType::gaussian);
                this->trainer->setCacheSize(100);
                this->trainer->setTolerance(0.000001);
                
                this->trainer->setProblem(prob);
                
                this->trainer->setNU((1 / (prob->l * C)));
                this->C = C;
                
                this->trainer->setKernelWidth(w);
                this->quadratic = NULL;
                this->clusterizeBSV = true;
        }
                
        SVClustering::~SVClustering() {
                delete trainer;
                delete eucDist;
                // if (quadratic != NULL) free(quadratic);
        }

        void SVClustering::learn() {
                free(quadratic);
                quadratic = NULL;
                this->numberOfClusters = -1;
                this->pointPerCluster.clear();
                this->labels.clear();
                this->trainer->learn();
        }

        void SVClustering::learn(DataSet *train) {
                free(quadratic);
                quadratic = NULL;
                this->numberOfClusters = -1;
                this->pointPerCluster.clear();
                this->labels.clear();
                this->trainer->learn(train);
        }
        
        void SVClustering::learn(struct svm_problem *p) {
                free(quadratic);
                quadratic = NULL;
                this->numberOfClusters = -1;
                this->pointPerCluster.clear();
                this->labels.clear();
                this->trainer->learn(p);
        }
        
                
        
        double SVClustering::getSphereRadius() {
                calculateSphereRadius();
                return sphereRadius;
        }
        
        void SVClustering::setTrainingSet(DataSet *t) {
                this->trainer->setTrainingSet(t);
        }
        
        DataSet *SVClustering::getTrainingSet() {
                return this->trainer->getTrainingSet();
        }
        
        
        void SVClustering::setProblem(struct svm_problem *p) {
                this->trainer->setProblem(p);
        }
        
        struct svm_problem *SVClustering::getProblem() {
                return this->trainer->getProblem();
        }
        
        
        struct svm_model *SVClustering::getModel() {
                return this->trainer->getModel();
        }

        struct svm_parameter *SVClustering::getParameters() {
                return trainer->getParameters();
        }
        
        double SVClustering::initialKernelWidth() {
                struct svm_problem *p = this->trainer->getProblem();
                int N = this->trainer->getProblem()->l;
                
                double max = eucDist->calculateDistance(p->x[0], p->x[1]);
                double dist;
                
                for (int i = 0; i < N; i++) {
                        for (int j = 0; j < N; j++) {
                                dist = eucDist->calculateDistance(p->x[i], p->x[j]);
                                if (dist > max) {
                                        max = dist;
                                }
                        }
                }
                
                this->trainer->setKernelWidth(1/max);
                
                return (1/max);
        }
        
        
        void SVClustering::setKernelWidth(double w) {
                this->trainer->setKernelWidth(w);
        }
        
        
        double SVClustering::getKernelWidth() {
                return this->trainer->getKernelWidth();
        }

        void SVClustering::setKernelType(int k) {
                this->trainer->setKernel(k);
        }       
        
        int SVClustering::getKernelType() {
                return this->trainer->getKernel();
        }
                
        void SVClustering::setSoftConstraint(double C) {
                this->C = C;
                
                /*      The OneClass SVM needs the NU param to be set in function of C param */
                this->trainer->setNU((1 / (this->trainer->getProblem()->l * C)));       
        }
        
        double SVClustering::getSoftConstraint() {
                return this->C;
        }


        double SVClustering::getSoftConstraintEstimate() {
                return C_estimation;
        }
        
        
        void SVClustering::clusterize() {
                this->separateClusters();
        }
        
        std::vector<int>& SVClustering::getClustersAssignment() {
                return labels;
        }
        
        std::vector<unsigned long int>& SVClustering::getPointPerCluster() {
                if (numberOfClusters > -1) {
                        return pointPerCluster;
                }
                
                // calculating points per cluster and number of clusters detected
                
                // theoretically impossible
                // here to detect some programming errors
                if (labels.size() == 0) {
                        this->numberOfClusters = 0;
                        return pointPerCluster;
                }
                
                std::vector<unsigned long int>::size_type idx;
                this->pointPerCluster.clear();
                this->pointPerCluster.push_back(0);
                this->numberOfClusters = 1;
        for (idx = 0; idx < labels.size(); idx++) {
                        if ((int)(pointPerCluster.size() - 1) < labels[idx]) {
                        this->numberOfClusters++;
                        this->pointPerCluster.push_back(1);
                }
                else {
                        this->pointPerCluster[labels[idx]]++;
                }
        }
                        
                return pointPerCluster;
        }
        
        unsigned long int SVClustering::getNumberOfClusters() {
                
                if (numberOfClusters > -1)
                        return numberOfClusters;
                
                // calculating points per cluster and number of clusters detected
                
                // theoretically impossible
                // here to detect some programming errors
                if (labels.size() == 0) {
                        this->numberOfClusters = 0;
                        return numberOfClusters;
                }
                
                
                std::vector<int>::size_type idx;
                this->pointPerCluster.push_back(0);
                this->numberOfClusters = 1;
        for (idx = 0; idx < labels.size(); idx++) {
                        if ((int)(pointPerCluster.size() - 1) < labels[idx]) {
                        this->numberOfClusters++;
                        this->pointPerCluster.push_back(1);
                }
                else {
                        this->pointPerCluster[labels[idx]]++;
                }
        }

                return numberOfClusters;
        }
        
        unsigned long int SVClustering::getNumberOfNonSingletonClusters() {
                return getNumberOfValidClusters(2);
        }

        unsigned long int SVClustering::getNumberOfValidClusters(unsigned long int threshold) {
                getPointPerCluster();
                unsigned long int counter = 0;
                for (unsigned long int i = 0; i < pointPerCluster.size(); i++) {
                        if (pointPerCluster[i] >= threshold) {
                                counter++;
                        }
                }
                return counter;
        }
                
        std::vector<int>& SVClustering::getOriginalClasses(int classOffset) {
                double *classes = this->trainer->getProblem()->y;
                unsigned long int i, N = this->trainer->getProblem()->l;
                
                originalClasses.resize(N);
                
                for (i = 0; i < N; i++) {
                        originalClasses[i] = (int)classes[i] - classOffset;
                }
                
                return originalClasses;
        }
        

        void SVClustering::setEuclideanDistance(EuclideanDistance *ed) {
                eucDist = ed;
        }       


        void SVClustering::useSpecialClusteringPolicyForBSV() {
                this->clusterizeBSV = true;
        }
        
        void SVClustering::useClassicClusteringPolicyForBSV() {
                this->clusterizeBSV = false;
        }

        bool SVClustering::usingSpecialClusteringPolicyForBSV() {
                return this->clusterizeBSV;
        }

        
        /*      PROTECTED METHODS FOLLOW        ************************************************************************/
        
        
         double SVClustering::beta(int i) {
                return this->trainer->getModel()->sv_coef[0][i] * C; 
         }
         
                 
        double SVClustering::beta(int i, bool scaled) {
                if (scaled) {
                        return this->trainer->getModel()->sv_coef[0][i];
                }
                else {
                        return this->trainer->getModel()->sv_coef[0][i] * getSoftConstraint();
                }
        }
        
         
         double SVClustering::rho() {
                return this->trainer->getModel()->rho[0] * C;
         }

         double SVClustering::rho(bool scaled) {
                if (scaled) {
                        return this->trainer->getModel()->rho[0];
                }
                else {
                        return this->trainer->getModel()->rho[0] * getSoftConstraint();
                }
         }
         
        
        void SVClustering::calculateQuadraticPartOfDistanceFromCenter() {
                
                if (quadratic != NULL) {
                        return;
                }
                
                struct svm_node **SV = this->trainer->getModel()->SV;
                int nSV = this->trainer->getModel()->l;
                double tmp = 0;
                
                
                for (int i = 0; i < nSV; i++) {
                        for (int j = 0; j < nSV; j++) {
                                tmp += beta(i) * beta(j) * kernelfunction(SV[i], SV[j], this->trainer->getModel()->param);
                        }
                }
                
                this->quadratic = (double *)malloc(sizeof(double));
                *(this->quadratic) = tmp;
        } 
         
        double SVClustering::calculateDistanceFromCenter(struct svm_node *x) {
                double sum = 0;
                
                struct svm_node **SV = this->trainer->getModel()->SV;
                int nSV = this->trainer->getModel()->l;
                
                
                for (int i = 0; i < nSV; i++) {
                        sum += beta(i) * kernelfunction(SV[i], x, this->trainer->getModel()->param);
                }
                
                
                calculateQuadraticPartOfDistanceFromCenter();
                
                // 1 is because (in case of RBF kernel used here) K(x,x) is equal to 1 
                // (therotically provable and pratically tested with libsvm) 
                //return sqrt(1 - (2 * sum) + *(this->quadratic));
                return 1 - (2 * sum) + *(this->quadratic);
        }
        
        void SVClustering::calculateSphereRadius() {
                
                // first tecnique
                //
                // assuming SVs have all the same distanace from center (theoretically a correct assumption)
                // to avoid a check on the distance of all SVs
                //
                this->sphereRadius = this->calculateDistanceFromCenter(this->trainer->getModel()->SV[0]);
                
                //if (this->sphereRadius >= 1) {
                //      free(quadratic);
                //      this->sphereRadius = this->calculateDistanceFromCenter(this->trainer->getModel()->SV[0]);
                //}

                // second tecnique
                //
                // we can use the average among SVs' distances to calculate the radius
                
                //int nSV = this->trainer->getModel()->l;
                //this->sphereRadius = 0;
                //for (int i = 0; i < nSV; i++) {
                //      this->sphereRadius += this->calculateDistanceFromCenter(this->trainer->getModel()->SV[i]);                      
                //}
                //this->sphereRadius = this->sphereRadius / nSV;

                //      third tecnique
                //      Using One Class SVM, Radius = sqrt(1 - rho/2)
                //sphereRadius = sqrt(1 - (rho()/2));
        }

        
        struct svm_node *SVClustering::pointOnThePath(struct svm_node *start, struct svm_node *stop, double interval) {
                int start_dim = 0, stop_dim = 0;
                struct svm_node *t, *z;
                
                t = start;
                while (t->index != -1) {
                        start_dim++;
                        t++;
                }
                
                t = stop;
                while (t->index != -1) {
                        stop_dim++;
                        t++;
                }
                
                start_dim++;
                stop_dim++;
                
                
                if (start_dim >= stop_dim) {
                        z = (struct svm_node *) malloc(sizeof(struct svm_node) * start_dim);
                }
                else {
                        z = (struct svm_node *) malloc(sizeof(struct svm_node) * stop_dim);
                }

                t = z;
                
                struct svm_node *x = start;
                struct svm_node *y = stop;
                                
                while(x->index != -1 && y->index !=-1) {
                        if (x->index == y->index) {
                                z->index = x->index;
                                z->value = ((y->value - x->value) * interval) + x->value;
                                x++;
                                y++;
                                z++;
                        }
                        else {
                                if (x->index > y->index) {      
                                        z->index = y->index;
                                        z->value = y->value * interval; // (y->value - 0) * interval + 0 
                                        ++y;
                                        ++z;
                                }
                                else {
                                        z->index = x->index;
                                        z->value = -(x->value * interval) + x->value;  // (0 - x->value) * interval + x->value 
                                        ++x;
                                        ++z;
                                }
                        }
                }

                while(x->index != -1) {
                        z->index = x->index;
                        z->value = -(x->value * interval) + x->value;  // (0 - x->value) * interval + x->value
                        ++x;
                        ++z;
                }

                while(y->index != -1) {
                        z->index = y->index;
                        z->value = y->value * interval; // (y->value - 0) * interval + 0 
                        ++y;
                        ++z;
                }
                
                z->index = -1;
                z->value = 0;
                
                return t;
        }
        
        
        void SVClustering::separateClusters() {
                static struct svm_node *z;
                bool flag;
                int i,j;
                double zdist;
                int N = this->trainer->getProblem()->l;
                
                std::pair<UndirectedGraph::edge_descriptor, bool> temp;
                        
                UndirectedGraph aGraph(N);
                
                // calculating adj matrix
                
                for (i = 0; i < N; i++) {
                        for (j = 0; j < N; j++) {
                                
                                // if the j is adjacent to i, then all j adjacent's are also adjacent to i.
                                if (j < i) {
                                        temp = edge(i, j, aGraph);
                                        if (temp.second == true) {
                                                for (int k = 0; k < N; k++) {
                                                        if (edge(i, k, aGraph).second || edge(j, k, aGraph).second) {
                                                                add_edge(i, k, aGraph);
                                                        }
                                                        else {
                                                                remove_edge(i, k, aGraph);
                                                        }
                                                }
                                        }               
                                }
                                else {
                                        // the i-th node is adiacent to itself - trivial
                                        if (i == j) {
                                                add_edge(i, j, aGraph);
                                                continue;
                                        }

                                        // we avoid to compute the same edge more than once                                     
                                        temp = edge(i, j, aGraph);
                                        if (temp.second == false) {
                                                flag = true;

                                                //      sampling the point on the path between x[i] and x[j]
                                                //      usually 10-20 points are sampled in order to limit running time
                                                //      here 10 points are sampled
                                                for (double k = 1; k <= 10; k = k + 1) {
                                                        z = pointOnThePath(this->trainer->getProblem()->x[i], this->trainer->getProblem()->x[j], k / 10);
                                                        
                                                        zdist = calculateDistanceFromCenter(z);

                                                        if (zdist > this->getSphereRadius()) {
                                                                flag = false;
                                                                break;
                                                        }
                                                }
                                                
                                                if (flag) {
                                                        add_edge(i, j, aGraph);
                                                        add_edge(j, i, aGraph);
                                                }
                                                else {
                                                        remove_edge(i, j, aGraph);
                                                        remove_edge(j, i, aGraph);
                                                }
                                        }
                                }
                        }
                }
                        
                // finding connected components (and so the clusters)
                
                labels.resize(N);
                connected_components(aGraph, &labels[0]);
                
                
                // clusterize the BSVs
                
                int nBSV = this->trainer->getModel()->lbounded;
                int *BSV_index = this->trainer->getModel()->BSV_index;

                int minimizer;
                double min_BSV_distance, d;
                
                for (int i = 0; i < nBSV; i++) {
                        min_BSV_distance = MAXFLOAT; 
                        minimizer = 0;
                        for (int j = 0; j < (int) labels.size(); j++) {
                                if (j == BSV_index[i]) continue;        //      ignore if the point is the current BSV
                                if (beta(j, true) == 1) continue;       //      ignore all BSVs
                                
                                d = eucDist->calculateDistance(this->trainer->getProblem()->x[BSV_index[i]], this->trainer->getProblem()->x[j]);
                                if (d < min_BSV_distance) {
                                        min_BSV_distance = d;
                                        minimizer = i;
                                }
                        }
                        labels[BSV_index[i]] = labels[minimizer];
                }
                
        }
        
}

---