KNN.cpp

Go to the documentation of this file.

/*
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or
 * (at your option) any later version.
 *
 * Written (W) 2006 Christian Gehl
 * Written (W) 2006-2009 Soeren Sonnenburg
 * Written (W) 2011 Sergey Lisitsyn
 * Written (W) 2012 Fernando José Iglesias García, cover tree support
 * Copyright (C) 2011 Berlin Institute of Technology and Max-Planck-Society
 */

#include <shogun/multiclass/KNN.h>
#include <shogun/labels/Labels.h>
#include <shogun/labels/MulticlassLabels.h>
#include <shogun/mathematics/Math.h>
#include <shogun/lib/Signal.h>
#include <shogun/lib/JLCoverTree.h>
#include <shogun/lib/Time.h>
#include <shogun/base/Parameter.h>
#include <shogun/multiclass/tree/KDTree.h>
#include <shogun/mathematics/eigen3.h>

#ifdef HAVE_CXX11
#include <shogun/lib/external/falconn/lsh_nn_table.h>
#endif

//#define DEBUG_KNN

using namespace shogun;
using namespace Eigen;

CKNN::CKNN()
: CDistanceMachine()
{
    init();
}

CKNN::CKNN(int32_t k, CDistance* d, CLabels* trainlab, KNN_SOLVER knn_solver)
: CDistanceMachine()
{
    init();

    m_k=k;

    ASSERT(d)
    ASSERT(trainlab)

    set_distance(d);
    set_labels(trainlab);
    m_train_labels.vlen=trainlab->get_num_labels();
    m_knn_solver=knn_solver;
}

void CKNN::init()
{
    /* do not store model features by default (CDistanceMachine::apply(...) is
     * overwritten */
    set_store_model_features(false);

    m_k=3;
    m_q=1.0;
    m_num_classes=0;
    m_leaf_size=1;
    m_knn_solver=KNN_BRUTE;
#ifdef HAVE_CXX11
    m_lsh_l = 0;
    m_lsh_t = 0;
#endif

    /* use the method classify_multiply_k to experiment with different values
     * of k */
    SG_ADD(&m_k, "m_k", "Parameter k", MS_NOT_AVAILABLE);
    SG_ADD(&m_q, "m_q", "Parameter q", MS_AVAILABLE);
    SG_ADD(&m_num_classes, "m_num_classes", "Number of classes", MS_NOT_AVAILABLE);
    SG_ADD(&m_leaf_size, "m_leaf_size", "Leaf size for KDTree", MS_NOT_AVAILABLE);
    SG_ADD((machine_int_t*) &m_knn_solver, "m_knn_solver", "Algorithm to solve knn", MS_NOT_AVAILABLE);
}

CKNN::~CKNN()
{
}

bool CKNN::train_machine(CFeatures* data)
{
    ASSERT(m_labels)
    ASSERT(distance)

    if (data)
    {
        if (m_labels->get_num_labels() != data->get_num_vectors())
            SG_ERROR("Number of training vectors does not match number of labels\n")
        distance->init(data, data);
    }

    SGVector<int32_t> lab=((CMulticlassLabels*) m_labels)->get_int_labels();
    m_train_labels=lab.clone();
    ASSERT(m_train_labels.vlen>0)

    int32_t max_class=m_train_labels[0];
    int32_t min_class=m_train_labels[0];

    for (int32_t i=1; i<m_train_labels.vlen; i++)
    {
        max_class=CMath::max(max_class, m_train_labels[i]);
        min_class=CMath::min(min_class, m_train_labels[i]);
    }

    for (int32_t i=0; i<m_train_labels.vlen; i++)
        m_train_labels[i]-=min_class;

    m_min_label=min_class;
    m_num_classes=max_class-min_class+1;

    SG_INFO("m_num_classes: %d (%+d to %+d) num_train: %d\n", m_num_classes,
            min_class, max_class, m_train_labels.vlen);

    return true;
}

SGMatrix<index_t> CKNN::nearest_neighbors()
{
    //number of examples to which kNN is applied
    int32_t n=distance->get_num_vec_rhs();
    //distances to train data
    float64_t* dists=SG_MALLOC(float64_t, m_train_labels.vlen);
    //indices to train data
    index_t* train_idxs=SG_MALLOC(index_t, m_train_labels.vlen);
    //pre-allocation of the nearest neighbors
    SGMatrix<index_t> NN(m_k, n);

    distance->precompute_lhs();

    //for each test example
    for (int32_t i=0; i<n && (!CSignal::cancel_computations()); i++)
    {
        SG_PROGRESS(i, 0, n)

        //lhs idx 0..num train examples-1 (i.e., all train examples) and rhs idx i
        distances_lhs(dists,0,m_train_labels.vlen-1,i);

        //fill in an array with 0..num train examples-1
        for (int32_t j=0; j<m_train_labels.vlen; j++)
            train_idxs[j]=j;

        //sort the distance vector between test example i and all train examples
        CMath::qsort_index(dists, train_idxs, m_train_labels.vlen);

#ifdef DEBUG_KNN
        SG_PRINT("\nQuick sort query %d\n", i)
        for (int32_t j=0; j<m_k; j++)
            SG_PRINT("%d ", train_idxs[j])
        SG_PRINT("\n")
#endif

        //fill in the output the indices of the nearest neighbors
        for (int32_t j=0; j<m_k; j++)
            NN(j,i) = train_idxs[j];
    }

    distance->reset_precompute();

    SG_FREE(train_idxs);
    SG_FREE(dists);

    return NN;
}

CMulticlassLabels* CKNN::apply_multiclass(CFeatures* data)
{
    if (data)
        init_distance(data);

    //redirecting to fast (without sorting) classify if k==1
    if (m_k == 1)
        return classify_NN();

    ASSERT(m_num_classes>0)
    ASSERT(distance)
    ASSERT(distance->get_num_vec_rhs())

    int32_t num_lab=distance->get_num_vec_rhs();
    ASSERT(m_k<=distance->get_num_vec_lhs())

    CMulticlassLabels* output=new CMulticlassLabels(num_lab);

    //labels of the k nearest neighbors
    int32_t* train_lab=SG_MALLOC(int32_t, m_k);

    SG_INFO("%d test examples\n", num_lab)
    CSignal::clear_cancel();

    //histogram of classes and returned output
    float64_t* classes=SG_MALLOC(float64_t, m_num_classes);

    switch (m_knn_solver)
    {
    case KNN_BRUTE:
    {
        //get the k nearest neighbors of each example
        SGMatrix<index_t> NN = nearest_neighbors();

        //from the indices to the nearest neighbors, compute the class labels
        for (int32_t i=0; i<num_lab && (!CSignal::cancel_computations()); i++)
        {
            //write the labels of the k nearest neighbors from theirs indices
            for (int32_t j=0; j<m_k; j++)
                train_lab[j] = m_train_labels[ NN(j,i) ];

            //get the index of the 'nearest' class
            int32_t out_idx = choose_class(classes, train_lab);
            //write the label of 'nearest' in the output
            output->set_label(i, out_idx + m_min_label);
        }

        break;
    }
    case KNN_COVER_TREE: // Use cover tree
    {
        // m_q != 1.0 not supported with cover tree because the neighbors
        // are not retrieved in increasing order of distance to the query
        float64_t old_q = m_q;
        if ( old_q != 1.0 )
            SG_INFO("q != 1.0 not supported with cover tree, using q = 1\n")

        // From the sets of features (lhs and rhs) stored in distance,
        // build arrays of cover tree points
        v_array< CJLCoverTreePoint > set_of_points  =
            parse_points(distance, FC_LHS);
        v_array< CJLCoverTreePoint > set_of_queries =
            parse_points(distance, FC_RHS);

        // Build the cover trees, one for the test vectors (rhs features)
        // and another for the training vectors (lhs features)
        CFeatures* r = distance->replace_rhs( distance->get_lhs() );
        node< CJLCoverTreePoint > top = batch_create(set_of_points);
        CFeatures* l = distance->replace_lhs(r);
        distance->replace_rhs(r);
        node< CJLCoverTreePoint > top_query = batch_create(set_of_queries);

        // Get the k nearest neighbors to all the test vectors (batch method)
        distance->replace_lhs(l);
        v_array< v_array< CJLCoverTreePoint > > res;
        k_nearest_neighbor(top, top_query, res, m_k);

#ifdef DEBUG_KNN
        SG_PRINT("\nJL Results:\n")
        for ( int32_t i = 0 ; i < res.index ; ++i )
        {
            for ( int32_t j = 0 ; j < res[i].index ; ++j )
            {
                printf("%d ", res[i][j].m_index);
            }
            printf("\n");
        }
        SG_PRINT("\n")
#endif

        for ( int32_t i = 0 ; i < res.index ; ++i )
        {
            // Translate from indices to labels of the nearest neighbors
            for ( int32_t j = 0; j < m_k; ++j )
                // The first index in res[i] points to the test vector
                train_lab[j] = m_train_labels.vector[ res[i][j+1].m_index ];

            // Get the index of the 'nearest' class
            int32_t out_idx = choose_class(classes, train_lab);
            output->set_label(res[i][0].m_index, out_idx+m_min_label);
        }

        m_q = old_q;

        break;
    }
    case KNN_KDTREE:
    {
        CFeatures* lhs = distance->get_lhs();
        CKDTree* kd_tree = new CKDTree(m_leaf_size);
        kd_tree->build_tree(dynamic_cast<CDenseFeatures<float64_t>*>(lhs));
        SG_UNREF(lhs);

        CFeatures* query = distance->get_rhs();
        kd_tree->query_knn(dynamic_cast<CDenseFeatures<float64_t>*>(query), m_k);
        SGMatrix<index_t> NN = kd_tree->get_knn_indices();
        for (int32_t i=0; i<num_lab && (!CSignal::cancel_computations()); i++)
        {
            //write the labels of the k nearest neighbors from theirs indices
            for (int32_t j=0; j<m_k; j++)
                train_lab[j] = m_train_labels[ NN(j,i) ];

            //get the index of the 'nearest' class
            int32_t out_idx = choose_class(classes, train_lab);
            //write the label of 'nearest' in the output
            output->set_label(i, out_idx + m_min_label);
        }
        SG_UNREF(query);
        break;
    }
#ifdef HAVE_CXX11
    case KNN_LSH:
    {
        CDenseFeatures<float64_t>* features = dynamic_cast<CDenseFeatures<float64_t>*>(distance->get_lhs());
        std::vector<falconn::DenseVector<double>> feats;
        for(int32_t i=0; i < features->get_num_vectors(); i++)
        {
            int32_t len;
            bool free;
            float64_t* vec = features->get_feature_vector(i, len, free);
            falconn::DenseVector<double> temp = Map<VectorXd> (vec, len);
            feats.push_back(temp);
        }

        falconn::LSHConstructionParameters params 
            = falconn::get_default_parameters<falconn::DenseVector<double>>(features->get_num_vectors(),
                                   features->get_num_features(),
                                   falconn::DistanceFunction::EuclideanSquared,
                                   true);
        SG_UNREF(features);
        if (m_lsh_l && m_lsh_t)
            params.l = m_lsh_l;

        auto lsh_table = falconn::construct_table<falconn::DenseVector<double>>(feats, params);
        if (m_lsh_t)
            lsh_table->set_num_probes(m_lsh_t);

        CDenseFeatures<float64_t>* query_features = dynamic_cast<CDenseFeatures<float64_t>*>(distance->get_rhs());
        std::vector<falconn::DenseVector<double>> query_feats;

        SGMatrix<index_t> NN (m_k, query_features->get_num_vectors());
        for(int32_t i=0; i < query_features->get_num_vectors(); i++)
        {
            int32_t len;
            bool free;
            float64_t* vec = query_features->get_feature_vector(i, len, free);
            falconn::DenseVector<double> temp = Map<VectorXd> (vec, len);
            auto indices = new std::vector<int32_t> ();
            lsh_table->find_k_nearest_neighbors(temp, (int_fast64_t)m_k, indices);
            memcpy(NN.get_column_vector(i), indices->data(), sizeof(int32_t)*m_k);
            delete indices;
        }
        
        for (int32_t i=0; i<num_lab && (!CSignal::cancel_computations()); i++)
        {
            //write the labels of the k nearest neighbors from theirs indices
            for (int32_t j=0; j<m_k; j++)
                train_lab[j] = m_train_labels[ NN(j,i) ];

            //get the index of the 'nearest' class
            int32_t out_idx = choose_class(classes, train_lab);
            //write the label of 'nearest' in the output
            output->set_label(i, out_idx + m_min_label);
        }
        SG_UNREF(query_features);
        break;
    }
#endif /* HAVE_CXX11 */
    }

    SG_FREE(classes);
    SG_FREE(train_lab);

    return output;
}

CMulticlassLabels* CKNN::classify_NN()
{
    ASSERT(distance)
    ASSERT(m_num_classes>0)

    int32_t num_lab = distance->get_num_vec_rhs();
    ASSERT(num_lab)

    CMulticlassLabels* output = new CMulticlassLabels(num_lab);
    float64_t* distances = SG_MALLOC(float64_t, m_train_labels.vlen);

    SG_INFO("%d test examples\n", num_lab)
    CSignal::clear_cancel();

    distance->precompute_lhs();

    // for each test example
    for (int32_t i=0; i<num_lab && (!CSignal::cancel_computations()); i++)
    {
        SG_PROGRESS(i,0,num_lab)

        // get distances from i-th test example to 0..num_m_train_labels-1 train examples
        distances_lhs(distances,0,m_train_labels.vlen-1,i);
        int32_t j;

        // assuming 0th train examples as nearest to i-th test example
        int32_t out_idx = 0;
        float64_t min_dist = distances[0];

        // searching for nearest neighbor by comparing distances
        for (j=0; j<m_train_labels.vlen; j++)
        {
            if (distances[j]<min_dist)
            {
                min_dist = distances[j];
                out_idx = j;
            }
        }

        // label i-th test example with label of nearest neighbor with out_idx index
        output->set_label(i,m_train_labels.vector[out_idx]+m_min_label);
    }

    distance->reset_precompute();

    SG_FREE(distances);
    return output;
}

SGMatrix<int32_t> CKNN::classify_for_multiple_k()
{
    ASSERT(m_num_classes>0)
    ASSERT(distance)
    ASSERT(distance->get_num_vec_rhs())

    int32_t num_lab=distance->get_num_vec_rhs();
    ASSERT(m_k<=num_lab)

    int32_t* output=SG_MALLOC(int32_t, m_k*num_lab);

    //working buffer of m_train_labels
    int32_t* train_lab=SG_MALLOC(int32_t, m_k);

    //histogram of classes and returned output
    int32_t* classes=SG_MALLOC(int32_t, m_num_classes);

    SG_INFO("%d test examples\n", num_lab)
    CSignal::clear_cancel();
    
    switch (m_knn_solver)
    {
    case KNN_COVER_TREE: // Use cover tree
    {
        //allocation for distances to nearest neighbors
        float64_t* dists=SG_MALLOC(float64_t, m_k);

        // From the sets of features (lhs and rhs) stored in distance,
        // build arrays of cover tree points
        v_array< CJLCoverTreePoint > set_of_points  =
            parse_points(distance, FC_LHS);
        v_array< CJLCoverTreePoint > set_of_queries =
            parse_points(distance, FC_RHS);

        // Build the cover trees, one for the test vectors (rhs features)
        // and another for the training vectors (lhs features)
        CFeatures* r = distance->replace_rhs( distance->get_lhs() );
        node< CJLCoverTreePoint > top = batch_create(set_of_points);
        CFeatures* l = distance->replace_lhs(r);
        distance->replace_rhs(r);
        node< CJLCoverTreePoint > top_query = batch_create(set_of_queries);

        // Get the k nearest neighbors to all the test vectors (batch method)
        distance->replace_lhs(l);
        v_array< v_array< CJLCoverTreePoint > > res;
        k_nearest_neighbor(top, top_query, res, m_k);

        for ( int32_t i = 0 ; i < res.index ; ++i )
        {
            // Handle the fact that cover tree doesn't return neighbors
            // ordered by distance

            for ( int32_t j = 0 ; j < m_k ; ++j )
            {
                // The first index in res[i] points to the test vector
                dists[j]     = distance->distance(res[i][j+1].m_index,
                            res[i][0].m_index);
                train_lab[j] = m_train_labels.vector[
                            res[i][j+1].m_index ];
            }

            // Now we get the indices to the neighbors sorted by distance
            CMath::qsort_index(dists, train_lab, m_k);

            choose_class_for_multiple_k(output+res[i][0].m_index, classes,
                    train_lab, num_lab);
        }

        SG_FREE(dists);
        break;
    }
    case KNN_KDTREE:
    {
        //allocation for distances to nearest neighbors
        float64_t* dists=SG_MALLOC(float64_t, m_k);

        CFeatures* lhs = distance->get_lhs();
        CKDTree* kd_tree = new CKDTree(m_leaf_size);
        kd_tree->build_tree(dynamic_cast<CDenseFeatures<float64_t>*>(lhs));
        SG_UNREF(lhs);

        CFeatures* data = distance->get_rhs();
        kd_tree->query_knn(dynamic_cast<CDenseFeatures<float64_t>*>(data), m_k);
        SGMatrix<index_t> NN = kd_tree->get_knn_indices();
        for (int32_t i=0; i<num_lab && (!CSignal::cancel_computations()); i++)
        {
            //write the labels of the k nearest neighbors from theirs indices
            for (int32_t j=0; j<m_k; j++)
            {
                train_lab[j] = m_train_labels[ NN(j,i) ];
                dists[j] = distance->distance(i, NN(j,i));
            }
            CMath::qsort_index(dists, train_lab, m_k);

            choose_class_for_multiple_k(output+i, classes, train_lab, num_lab);
        }
        break;
    }
    default:
    {
        //get the k nearest neighbors of each example
        SGMatrix<index_t> NN = nearest_neighbors();

        for (int32_t i=0; i<num_lab && (!CSignal::cancel_computations()); i++)
        {
            //write the labels of the k nearest neighbors from theirs indices
            for (int32_t j=0; j<m_k; j++)
                train_lab[j] = m_train_labels[ NN(j,i) ];

            choose_class_for_multiple_k(output+i, classes, train_lab, num_lab);
        }

    }

    }

    SG_FREE(train_lab);
    SG_FREE(classes);

    return SGMatrix<int32_t>(output,num_lab,m_k,true);
}

void CKNN::init_distance(CFeatures* data)
{
    if (!distance)
        SG_ERROR("No distance assigned!\n")
    CFeatures* lhs=distance->get_lhs();
    if (!lhs || !lhs->get_num_vectors())
    {
        SG_UNREF(lhs);
        SG_ERROR("No vectors on left hand side\n")
    }
    distance->init(lhs, data);
    SG_UNREF(lhs);
}

bool CKNN::load(FILE* srcfile)
{
    SG_SET_LOCALE_C;
    SG_RESET_LOCALE;
    return false;
}

bool CKNN::save(FILE* dstfile)
{
    SG_SET_LOCALE_C;
    SG_RESET_LOCALE;
    return false;
}

void CKNN::store_model_features()
{
    CFeatures* d_lhs=distance->get_lhs();
    CFeatures* d_rhs=distance->get_rhs();

    /* copy lhs of underlying distance */
    distance->init(d_lhs->duplicate(), d_rhs);

    SG_UNREF(d_lhs);
    SG_UNREF(d_rhs);
}

int32_t CKNN::choose_class(float64_t* classes, int32_t* train_lab)
{
    memset(classes, 0, sizeof(float64_t)*m_num_classes);

    float64_t multiplier = m_q;
    for (int32_t j=0; j<m_k; j++)
    {
        classes[train_lab[j]]+= multiplier;
        multiplier*= multiplier;
    }

    //choose the class that got 'outputted' most often
    int32_t out_idx=0;
    float64_t out_max=0;

    for (int32_t j=0; j<m_num_classes; j++)
    {
        if (out_max< classes[j])
        {
            out_idx= j;
            out_max= classes[j];
        }
    }

    return out_idx;
}

void CKNN::choose_class_for_multiple_k(int32_t* output, int32_t* classes, int32_t* train_lab, int32_t step)
{
    //compute histogram of class outputs of the first k nearest neighbours
    memset(classes, 0, sizeof(int32_t)*m_num_classes);

    for (int32_t j=0; j<m_k; j++)
    {
        classes[train_lab[j]]++;

        //choose the class that got 'outputted' most often
        int32_t out_idx=0;
        int32_t out_max=0;

        for (int32_t c=0; c<m_num_classes; c++)
        {
            if (out_max< classes[c])
            {
                out_idx= c;
                out_max= classes[c];
            }
        }

        output[j*step]=out_idx+m_min_label;
    }
}