svc.c

Go to the documentation of this file.

#include <getopt.h>     //      command line parsing
#include <fcntl.h>      //      file handling

#include <errno.h>      //      standard C error handling
#include <err.h>        //      standard C error handling
#include <cstdio>       //      standard C I/O functions
#include <cstdlib>      //      standard C library functions
#include <cstring>      //      standard C string
#include <cmath>        //      standard C math functions

#include <iostream>     //      C++ I/O functions
#include <fstream>      //      C++ I/O functions
#include <string>       //      C++ string

#include <boost/timer.hpp>      //      Boost Library - Timer

#include "svm.h"                                //      SVM Library
#include "LibsvmFileLoader.h"   //      File Loader
#include "SVClustering.h"               //      Classic Support Vector Clustering with Complete Graph Cluster Labeling
#include "CCLSVClustering.h"    //      SVC with Cone Cluster Labeling
#include "GKWGenerator.h"               //      Q Generator
#include "KernelMachines.h"             //      Kernel Machine utilities


#include "L1Distance.h"                 //      L1 distance 
#include "L2Distance.h"                 //      L2 distance


#define DEFAULT_KW                              -1              //      default kernel width (no specified - rely on Kernel Width Generator)
#define DEFAULT_C                               1               //      default soft constraint value
#define DEFAULT_ALG                             0               //      default cluster labeling algorithm
#define DEFAULT_RUNS                    100000  //      default number of runs
#define DEFAULT_L_OFFSET                1               //      default label offset
#define DEFAULT_VALID_THRESHOLD 2               //      default minimum number of elements per cluster
#define DEFAULT_SOFTENING               1               //      default softening constant for kernel width generation
#define DEFAULT_EDIST                   2               //      default euclidean distance value (no specified - rely on SVC constructor)
#define DEFAULT_KERNEL                  -1              //      default kernel type (no specified - rely on SVC constructor)

void printAllPoints(damina::SVClustering *, unsigned long int, char*);

void countRightAndWrongClassificationsPerClass(damina::SVClustering *, vector<int> &, vector<int> &, vector<int> &, vector<int> &, int, int, vector<unsigned long int> &, unsigned long int);

void computeEvaluatingMeasreusPerClass(vector<int>, vector<int>, vector<int>, vector<int>, vector<double> &, vector<double> &, vector<double> &, vector<double> &, vector<double> &, double &, double &);
                
struct svc_options *parseCommandLineArgs(int, char **);

void usage(char **);

void printUsedCommandLine(int, char *[]);

enum {CCL = 0, CG = 1}; /* alg_type: cluster labeling algorithms types */


enum {L1 = 1, L2 = 2};  

enum {BSV_CLASSIC = 0, BSV_ENHANCED = 1};


struct svc_options {
        unsigned long int runs;                         //      number of runs
        unsigned long int alg_type;                     //      algorithm type
        double  q;                                                      //      gaussian kernel width
        double C;                                                       //      soft margin constant
        char *filename;                                         //      data-set filename
        unsigned long int clusters;                     //      number of clusters
        unsigned long int l_offset;                     //      labels offset
        unsigned long int valid_threshold;      //      minimal clusters cardinality 
        double softening;                                       //      softening constant for kernel width generation
        char *detfile;                                          //      file to store the detailed cluster assignment
        unsigned long int e_dist;                       //      euclidean distance type
        int kernel;                                                     //      kernel type
        unsigned int bsv;                                               //      BSV clustering policy
};

int main(int argc, char *argv[]) {

        using namespace std;
        using namespace damina;
        using namespace boost;
        
        struct svc_options *opt;
        
        opt = parseCommandLineArgs(argc, argv); //      parse command line to get parameters
        if (opt == NULL) {                                              //      parsing fails
                usage(argv);                                            //      show usage
                exit(1);
        }

        LibsvmFileLoader *fl;
        string datafile(opt->filename);
        struct svm_problem *p = NULL;
        cout << "==================================" << endl;
        cout << "Loading file" << endl;
        cout << "==================================" << endl;
        fl = new LibsvmFileLoader(datafile);
        timer tfl;
        p = fl->load(false);
        cout << "File loaded in " << tfl.elapsed() << " secs" << endl;
        cout << "Number of objects: " << p->l << endl << endl;
        delete fl;

        
        SVClustering *svc;                      //      Support Vector Clustering Engine
        switch (opt->alg_type) {
                case CG: svc = new SVClustering(opt->C, p); break;              //      classic one
                case CCL: svc = new CCLSVClustering(opt->C, p); break;  //      with CCL
                default:
                        svc = new SVClustering(opt->C, p); break;               //      default is classic
        }

        if (opt->bsv) {
                svc->useSpecialClusteringPolicyForBSV();
        }
        else {
                svc->useClassicClusteringPolicyForBSV();
        }
                
        if (opt->kernel != DEFAULT_KERNEL) {
                svc->setKernelType(opt->kernel);
        }

        GKWGenerator *kwg;
        if (opt->q == -1) {     //      no kernel width value provided on command line  
                kwg = new GKWGenerator(svc, 0.000000001);
        }
        else {  //      kernel width value provided on command line
                kwg = new GKWGenerator(svc, 0.000000001, opt->q);
        }

        kwg->setGrowthSoftening(opt->softening);

        switch (opt->e_dist) {
                case L1: { svc->setEuclideanDistance(new L1Distance()); kwg->setEuclideanDistance(new L1Distance()); } break;
                case L2: { svc->setEuclideanDistance(new L2Distance()); kwg->setEuclideanDistance(new L2Distance()); } break;
                default: ;
        }

        bool target = false;
        unsigned long int invalid = 0, runs = 0;
        double q, prev_q = -1, prev_C = -1;
        
        unsigned long int nrClusters, 
                                  // nrNonSingletonClusters;
                                  nrValidClusters;
        
        timer partial,  //      timer for intermediate steps
                        total;  //      timer for overall time measurement 
        
        vector<unsigned long int> ppc;
        vector<int> tp, fp, tn, fn;
        vector<double> precision, recall, f1, contamination, specificity;
        double ma, accuracy, overall_time = 0;
        
        while (true) {
                prev_q = q;
                prev_C = svc->getSoftConstraintEstimate();
                
                runs++;
                if (runs > opt->runs) break;

                q = kwg->getNextKernelWidthValue();
                if (q == -1) {
                        cout << "No more valid kernel width values" << endl << endl;
                        break;
                }
                
                if (opt->runs > 1) {
                        cout << "==================================" << endl;
                        cout << " Finding Kernel Width # " << runs << endl;
                        cout << "==================================" << endl;
                }

                svc->setKernelWidth(q);
                
                total.restart();                                //      start total time counter
                partial.restart();                              //      start domain description process time counter

                svc->learn();                                   //      domain description process
                cout << "radius = " << svc->getSphereRadius() << endl;
                cout << endl << "Domain description time taken: ";
                cout << partial.elapsed() << " seconds " << endl;               //      stop domain description process time counter

                partial.restart();                              //      start clustering process time counter

                svc->clusterize();                              //      cluster labeling process
                cout << "Clustering time taken: ";
                cout << partial.elapsed() << " seconds " << endl;       //      stop clustering process time counter

                cout << "Total time taken: ";
                overall_time += total.elapsed();
                cout << total.elapsed() << " seconds " << endl << endl;         //      stop total time counter
                
                
                nrClusters = svc->getNumberOfClusters();
                // nrNonSingletonClusters = svc->getNumberOfNonSingletonClusters();
                nrValidClusters = svc->getNumberOfValidClusters(opt->valid_threshold);

                cout << "Soft Constraint Estimate: " << svc->getSoftConstraintEstimate() << endl;       
                cout << "Current Soft Constraint: " << svc->getSoftConstraint() << endl;
                cout << "Kernel Width: " << q << endl;
                cout << "Number of clusters: " << nrClusters << endl;
                cout << "Number of valid clusters: " << nrValidClusters << endl;
                cout << "Number of non valid clusters (cardinality < " << opt->valid_threshold << "): " << nrClusters - nrValidClusters << endl;
                
                cout << "Points per cluster: " << endl;
        
                ppc.clear();    
                ppc = svc->getPointPerCluster();
        
                for (unsigned long int i = 0; i < ppc.size(); i++) {
                        if (ppc[i] >= opt->valid_threshold) {
                                cout << "\tCluster " << i << ": " << ppc[i] << endl;
                        }
                }
                cout << endl << endl;

                if (opt->detfile != NULL) {
                        printAllPoints(svc, opt->l_offset, opt->detfile);
                }

                countRightAndWrongClassificationsPerClass(svc, tp, fp, tn, fn, opt->clusters, opt->l_offset, ppc, opt->valid_threshold);
                computeEvaluatingMeasreusPerClass(tp, fp, tn, fn, contamination, specificity, precision, recall, f1, ma, accuracy);
                
                for (unsigned long int i = 0; i < tp.size(); i++) {
                        cout << "Class " << i << endl << "\tTP: " << tp[i] << "\tFP: " << fp[i] << endl << "\tFN: " << fn[i] << "\tTN: " << tn[i] << endl << endl; 
                        cout << "Precision: " << precision[i] << " - Recall: " << recall[i] << " - F1: " << f1[i] << endl;
                        cout << "Completeness: " << recall[i] << " - Contamination: " << contamination[i] << endl;
                        cout << "Sensitivity: " << recall[i]  << " - Specificity: " << specificity[i] << endl << endl;
                }
                
                if (opt->clusters > 2) {
                        cout << "Macroaveraging: " << ma << endl; 
                }
                
                cout << "Accuracy: " << accuracy << endl << endl;

                // check stop criteria
                if (target) {
                        if (nrValidClusters != opt->clusters) {
                                break;
                        }
                        else {
                                if (nrClusters - nrValidClusters > invalid) {
                                        break;
                                }
                        }
                }
                else {
                        if (nrValidClusters == opt->clusters) {
                                target = true;
                                invalid = nrClusters - nrValidClusters;
                        }
                }
                        
        }

        cout << "==================================" << endl;
        cout << "Clustering process finished" << endl;
        cout << "==================================" << endl;
        if (opt->runs > 1) {
                cout << "Check out the last run and second last run." << endl;
                cout << "Kernel Width Found: " << prev_q << endl;
        }   
        cout << "Soft Constraint Estimate: " << prev_C << endl;
    cout << "Overall time elapsed: " << overall_time << " secs" << endl;
//      cout << "Memory used: " << sizeof(svc) + sizeof(kwg) << endl;
        cout << "Original command line: "; 
        printUsedCommandLine(argc, argv);

        delete svc;
        delete kwg;
        free(opt);
}


void printAllPoints(damina::SVClustering *clust, unsigned long int classOffset, char *filename) {
        
        filebuf fb;
        fb.open(filename, ios::app);
        ostream out(&fb);
                
        std::vector<int> labels = clust->getClustersAssignment();               //      clustering result
        std::vector<int> orig_c = clust->getOriginalClasses((int)classOffset);          //      original classification

        for (unsigned long int i = 0; i < labels.size(); i++) {
                out << "Point #" << i << " - " << "Cluster: " << labels[i] << ", Class: " << orig_c[i] << endl;
        }
        out << endl;
        fb.close();
}

void countRightAndWrongClassificationsPerClass(damina::SVClustering *clust, vector<int> &tp, vector<int> &fp, vector<int> &tn, vector<int> &fn, int nrClasses, int classOffset, vector<unsigned long int> &ppc, unsigned long int valid_t) {
        
        std::vector<int> labels = clust->getClustersAssignment();               //      clustering result
        std::vector<int> orig_c = clust->getOriginalClasses(classOffset);               //      original classification
        
        int datasetDimension = clust->getProblem()->l;

        tp.clear();     
        fp.clear();     
        tn.clear();     
        fn.clear();     
        
        tp.resize(nrClasses);
        fp.resize(nrClasses);
        tn.resize(nrClasses);
        fn.resize(nrClasses);   
        
        for (unsigned long int i = 0; i < labels.size(); i++) {
                if (orig_c[i] >= nrClasses) continue;
                
                if ((ppc[labels[i]] >= valid_t) && (labels[i] == orig_c[i])) {
                        tp[orig_c[i]]++;
                }
                else {
                        fp[labels[i]]++;
                        fn[orig_c[i]]++;
                }
        }
        
        for (int i = 0; i < nrClasses; i++) {
                tn[i] = datasetDimension - tp[i] - fp[i] - fn[i];
        }
}


void computeEvaluatingMeasreusPerClass(vector<int> tp, vector<int> fp, vector<int> tn, vector<int> fn, 
                        vector<double> &contamination, vector<double> &specificity, 
                        vector<double> &precision, vector<double> &recall, vector<double> &f1, double &ma, double &accuracy) {

        unsigned long int nrClasses = tp.size();
        int datasedDim = tp[0] + fp[0] + tn[0] + fn[0];
        
        precision.clear();
        recall.clear();
        f1.clear();
        
        contamination.resize(nrClasses);
        precision.resize(nrClasses);
        recall.resize(nrClasses);
        f1.resize(nrClasses);
        specificity.resize(nrClasses);
        
        ma = 0;
        accuracy = 0;
        
        for (unsigned long int i = 0; i < nrClasses ; i++) {
                accuracy += tp[i];
                precision[i] = ((double)tp[i] / (double)(tp[i] + fp[i])) * 100;
                recall[i] = ((double)tp[i] / (double)(tp[i] + fn[i])) * 100;
                specificity[i] = ((double)tn[i] / (double)(tn[i] + fp[i])) * 100;
                contamination[i] = (double)fp[i] / (double)tp[i];
                if (precision[i] + recall[i] == 0) {
                        f1[i] = -1; 
                }
                else {
                        f1[i] = (2 * (precision[i] * recall[i])) / (precision[i] + recall[i]);
                        ma += f1[i];
                }

        }
                
        ma = ma / (double) nrClasses;
        accuracy = accuracy / (double)datasedDim * 100;
}


struct svc_options *parseCommandLineArgs(int argc, char **argv) {

        using namespace damina;
        
        struct svc_options *opt; 
        opt = (struct svc_options *)malloc(sizeof(struct svc_options));
        
        opt->filename = (char *)malloc(2 * sizeof(char));
        strcpy(opt->filename, "");
        
        opt->q = DEFAULT_KW;
        opt->C = DEFAULT_C;
        opt->alg_type = DEFAULT_ALG;
        opt->runs = DEFAULT_RUNS;
        opt->l_offset = DEFAULT_L_OFFSET;
        opt->valid_threshold = DEFAULT_VALID_THRESHOLD;
        opt->softening = DEFAULT_SOFTENING;
        opt->detfile = NULL;
        opt->e_dist = DEFAULT_EDIST;
        opt->kernel = DEFAULT_KERNEL;
        opt->bsv = BSV_CLASSIC;

        bool file_f = true;
        bool cluster_f = false;
        
        int ch;
        
        /* options descriptor */
        static struct option longopts[] = {
                { "kernel-width",       required_argument,      NULL,   'q' },
                { "soft-margin",        required_argument,      NULL,   'C' },
                { "num-clusters",       required_argument,      NULL,   'c' },
                { "labeling-alg",       required_argument,      NULL,   'l' },
                { "dataset-file",       required_argument,      NULL,   'f' },
                { "runs",                       required_argument,      NULL,   'r' },
                { "labels-offset",      required_argument,      NULL,   'o' },
                { "valid-threshold",required_argument,  NULL,   't' },
                { "softening",          required_argument,      NULL,   's' },
                { "dump-details",       required_argument,      NULL,   'D' },
                { "euclidean-distance", required_argument,      NULL,   'e' },
                { "kernel-type",        required_argument,      NULL,   'k' },
                { "bsv-clustering",     required_argument,      NULL,   'b' }

        };

        while ((ch = getopt_long(argc, argv, "q:C:c:l:f:r:o:t:s:D:e:k:b:", longopts, NULL)) != -1) {
                switch (ch) {
                                
                        case 'b': {
                                        opt->bsv = (int)strtol(optarg, (char**)NULL, 10);
                                        //if ((errno == EINVAL) || (opt->bsv != (int)BSV_CLASSIC) || (opt->bsv != (int)BSV_ENHANCED))
                                        //      err(1, "%s is an invalid value for BSV Clustering Policy. It must be an inter value in [%d,%d].", optarg, BSV_CLASSIC, BSV_ENHANCED);
                                        break;
                        }


                        case 'k': {
                                        opt->kernel = (int)strtol(optarg, (char**)NULL, 10);
                                        //if ((errno == EINVAL) || (opt->kernel != KernelType::gaussian) || (opt->kernel != KernelType::laplacian))
                                        //      err(1, "%s is an invalid value for Kernel Type. It must be %d or %d.", optarg, KernelType::gaussian, KernelType::laplacian);
                                        break;
                        }

                        case 'e': {
                                        opt->e_dist = (int)strtol(optarg, (char**)NULL, 10);
                                        if ((errno == EINVAL) || (opt->e_dist < 1) || (opt->e_dist > 2))
                                                err(1, "%s is an invalid value for Euclidean Distance type. It must be an inter value in [1,2].", optarg);
                                        break;
                        }

                        case 'q': {
                                        opt->q = (double)strtod(optarg, (char**)NULL);
                                        if ((errno == EINVAL) || (opt->q <= 0))
                                                err(1, "%s is an invalid value for Kernel Width. It must be a real value greater than 0.", optarg);
                                        break;
                        }
                        
                        case 'C': {
                                        opt->C = (double)strtod(optarg, (char**)NULL);
                                        if ((errno == EINVAL) || (opt->C <= 0))
                                                err(1, "%s is an invalid value for Soft Margin Constant. It must be a real value greater than 0.", optarg);
                                        break;
                        }
                        
                        case 'c': {
                                        opt->clusters = (int)strtol(optarg, (char**)NULL, 10);
                                        if ((errno == EINVAL) || (opt->clusters <= 0))
                                                err(1, "%s is an invalid value for Number of Clusters. It must be an integer value greater than 0.", optarg);
                                        
                                        cluster_f = true;
                                        break;
                        }
                        
                        case 'r': {
                                        opt->runs = (int)strtol(optarg, (char**)NULL, 10);
                                        if ((errno == EINVAL) || (opt->runs <= 0))
                                                err(1, "%s is an invalid value for Number of Runs. It must be an integer value greater than 0.", optarg);
                                        break;
                        }
                        
                        case 's': {
                                        opt->softening = (double)strtod(optarg, (char**)NULL);
                                        if ((errno == EINVAL) || (opt->softening <= 0))
                                                err(1, "%s is an invalid value for the Softening Constant. It must be an real value greater than 0.", optarg);
                                        break;
                        }
                
                        case 't': {
                                        opt->valid_threshold = (int)strtol(optarg, (char**)NULL, 10);
                                        if ((errno == EINVAL) || (opt->valid_threshold <= 0))
                                                err(1, "%s is an invalid value for Minimum Clusters Cardinality. It must be an integer value greater than 0.", optarg);
                                        break;
                        }

                        case 'o': {
                                        opt->l_offset = (int)strtol(optarg, (char**)NULL, 10);
                                        if ((errno == EINVAL) || (opt->l_offset < 0))
                                                err(1, "%s is an invalid value for Labels Offset. It must be an integer value greater than or equal to 0.", optarg);
                                        break;
                        }
                        
                        case 'l': {
                                        opt->alg_type = (int)strtol(optarg, (char**)NULL, 10);
                                        if ((errno == EINVAL) || (opt->alg_type > 1) || (opt->alg_type < 0))
                                                err(1, "%s is an invalid value for Cluster Labeling algorithm. It must be 0 or 1.", optarg);
                                        break;
                        }
                        
                        case 'f': {
                                        if ((open(optarg, O_RDONLY, 0)) == -1) {
                        err(1, "Unable to open %s", optarg);
                        break;
                                        }                       

                                        free(opt->filename);    
                                        opt->filename = (char *)malloc((strlen(optarg) + 1) * sizeof(char));
                                        strcpy(opt->filename, optarg);
                                                                                        
                                        file_f = true;
                                        break;
                        }
                        
                        case 'D': {
                                        free(opt->detfile);
                                        opt->detfile = (char *)malloc((strlen(optarg) + 1) * sizeof(char));
                                        strcpy(opt->detfile, optarg);
                                        break;
                        }
                        
                        default: {
                                usage(argv);
                        }
                }
        }
        argc -= optind;
        argv += optind; 

        if ((file_f == false) || (cluster_f == false)) {
                return NULL;
        }
        return opt;
}

void usage(char **argv) {
        using namespace damina;
        
        cerr << endl << "Usage: " << argv[0] << " -q (0,+inf) -C (0,+inf) -b [" << BSV_CLASSIC << "," << BSV_ENHANCED << "] -k" << KernelType::gaussian << "," << KernelType::laplacian << " -e [1,2] -c [1, N] -l [0-1] -f <dataset-filename> -D <detail-filename> -r [1, N] -o [0, N] -t [1,N] -s (0,1]" << endl << endl;
        cerr << "\t -q: the kernel width value must be a real value greater than zero - DEFAULT: auto" << endl; 
        cerr << "\t -C: the soft margin constant must be a real value greater than zero - DEFAULT: " << DEFAULT_C << endl;
        cerr << "\t -b: the BSV Clustering Policy - DEFAULT: " << BSV_ENHANCED << endl;
                cerr << "\t\t " << BSV_CLASSIC << " for classic policy" << endl;
                cerr << "\t\t " << BSV_ENHANCED << " for enhanced policy (aim to get better results)" << endl;
        cerr << "\t -k: the Kernel Type - DEFAULT: " << KernelType::gaussian << endl;
                cerr << "\t\t " << KernelType::gaussian << " for Gaussian Kernel (K(x,y) = exp(-gamma * ||x-y||^2))" << endl;
                cerr << "\t\t " << KernelType::laplacian << " for Laplacian Kernel (K(x,y) = exp(-gamma*|x-y|_1)) " << endl;
        cerr << "\t -e: the Euclidean Measure to use in vector-distance measurement - DEFAULT: " << L2 << endl;
                cerr << "\t\t " << L1 << " for L1 distance (sum |x_i - y_i|) " << endl;
                cerr << "\t\t " << L2 << " for L2 distance (sum sqrt(||x_i - y_i||^2)) " << endl;
        cerr << "\t -l: the cluster labeling algorithm type - DEFAULT: " << DEFAULT_ALG << endl;
                cerr << "\t\t " << CCL << " for Cone Cluster Labeling" << endl;
                cerr << "\t\t " << CG << " for Complete Graph Cluster Labeling" << endl;
        cerr << "\t -c: the number of cluster to request must be an integer value greater than zero" << endl;
        cerr << "\t -f: the filename of the dataset file, in LIBSVM format" << endl;
        cerr << "\t -D: the filename of the file to store detailed clusters assignment" << endl;
        cerr << "\t -r: the number of runs to perform - DEFAULT: " << DEFAULT_RUNS      << endl;
        cerr << "\t -o: the labels offset (i.e. the label of the first class) - DEFAULT: " << DEFAULT_L_OFFSET << endl;
        cerr << "\t -s: the softening constant for kernel width generation - DEFAULT: " << DEFAULT_SOFTENING << endl;
        cerr << "\t -t: the minimum clusters cardinality - DEFAULT: " << DEFAULT_VALID_THRESHOLD << endl << endl;
}

void printUsedCommandLine(int argc, char *argv[]) {
        cout << argv[0] << " ";

        for (int i = 1; i <= argc; i++) {
                cout << argv[i] << " ";
        }
        cout << endl << endl;
}

---