Kernel.h

Go to the documentation of this file.

/*
 * EXCEPT FOR THE KERNEL CACHING FUNCTIONS WHICH ARE (W) THORSTEN JOACHIMS
 * COPYRIGHT (C) 1999  UNIVERSITAET DORTMUND - ALL RIGHTS RESERVED
 *
 * 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) 1999-2009 Soeren Sonnenburg
 * Written (W) 1999-2008 Gunnar Raetsch
 * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society
 */

#ifndef _KERNEL_H___
#define _KERNEL_H___

#include "lib/common.h"
#include "lib/Signal.h"
#include "lib/File.h"
#include "lib/Mathematics.h"
#include "base/SGObject.h"
#include "features/Features.h"
#include "kernel/KernelNormalizer.h"

#include <vector>
#include <set>
#include <string>

namespace shogun
{
    class CFile;
    class CFeatures;
    class CKernelNormalizer;
    enum EFeatureType;
    enum EFeatureClass;

#ifdef USE_SHORTREAL_KERNELCACHE
    typedef float32_t KERNELCACHE_ELEM;
#else
    typedef float64_t KERNELCACHE_ELEM;
#endif

typedef int64_t KERNELCACHE_IDX;


enum EOptimizationType
{
    FASTBUTMEMHUNGRY,
    SLOWBUTMEMEFFICIENT
};

enum EKernelType
{
    K_UNKNOWN = 0,
    K_LINEAR = 10,
    K_POLY = 20,
    K_GAUSSIAN = 30,
    K_GAUSSIANSHIFT = 32,
    K_GAUSSIANMATCH = 33,
    K_HISTOGRAM = 40,
    K_SALZBERG = 41,
    K_LOCALITYIMPROVED = 50,
    K_SIMPLELOCALITYIMPROVED = 60,
    K_FIXEDDEGREE = 70,
    K_WEIGHTEDDEGREE =    80,
    K_WEIGHTEDDEGREEPOS = 81,
    K_WEIGHTEDDEGREERBF = 82,
    K_WEIGHTEDCOMMWORDSTRING = 90,
    K_POLYMATCH = 100,
    K_ALIGNMENT = 110,
    K_COMMWORDSTRING = 120,
    K_COMMULONGSTRING = 121,
    K_SPECTRUMMISMATCHRBF = 122,
    K_COMBINED = 140,
    K_AUC = 150,
    K_CUSTOM = 160,
    K_SIGMOID = 170,
    K_CHI2 = 180,
    K_DIAG = 190,
    K_CONST = 200,
    K_DISTANCE = 220,
    K_LOCALALIGNMENT = 230,
    K_PYRAMIDCHI2 = 240,
    K_OLIGO = 250,
    K_MATCHWORD = 260,
    K_TPPK = 270,
    K_REGULATORYMODULES = 280,
    K_SPARSESPATIALSAMPLE = 290,
    K_HISTOGRAMINTERSECTION = 300
};

enum EKernelProperty
{
    KP_NONE = 0,
    KP_LINADD = 1,  // Kernels that can be optimized via doing normal updates w + dw
    KP_KERNCOMBINATION = 2, // Kernels that are infact a linear combination of subkernels K=\sum_i b_i*K_i
    KP_BATCHEVALUATION = 4  // Kernels that can on the fly generate normals in linadd and more quickly/memory efficient process batches instead of single examples
};

template <class T> struct K_THREAD_PARAM
{
    CKernel* kernel;
    int32_t start;
    int32_t end;
    int32_t total_start;
    int32_t total_end;
    int32_t m;
    int32_t n;
    T* result;
    bool symmetric;
    bool verbose;
};

class CSVM;

class CKernel : public CSGObject
{
    friend class CVarianceKernelNormalizer;
    friend class CSqrtDiagKernelNormalizer;
    friend class CAvgDiagKernelNormalizer;
    friend class CRidgeKernelNormalizer;
    friend class CFirstElementKernelNormalizer;
    friend class CMultitaskKernelNormalizer;
    friend class CMultitaskKernelMklNormalizer;
    friend class CMultitaskKernelMaskNormalizer;
    friend class CMultitaskKernelMaskPairNormalizer;
    friend class CTanimotoKernelNormalizer;
    friend class CDiceKernelNormalizer;
    friend class CZeroMeanCenterKernelNormalizer;

    public:

        CKernel();


        CKernel(int32_t size);

        CKernel(CFeatures* l, CFeatures* r, int32_t size);

        virtual ~CKernel();

        inline float64_t kernel(int32_t idx_a, int32_t idx_b)
        {
            if (idx_a<0 || idx_b<0 || idx_a>=num_lhs || idx_b>=num_rhs)
            {
                SG_ERROR("Index out of Range: idx_a=%d/%d idx_b=%d/%d\n",
                        idx_a,num_lhs, idx_b,num_rhs);
            }

            return normalizer->normalize(compute(idx_a, idx_b), idx_a, idx_b);
        }

        void get_kernel_matrix(float64_t** dst, int32_t* m, int32_t* n);


        virtual std::vector<float64_t> get_kernel_col(int32_t j)
        {

            std::vector<float64_t> col = std::vector<float64_t>(num_rhs);

            for (int32_t i=0; i!=num_rhs; i++)
            {
                col[i] = kernel(i,j);
            }

            return col;

        }


        virtual std::vector<float64_t> get_kernel_row(int32_t i)
        {

            std::vector<float64_t> row = std::vector<float64_t>(num_lhs);

            for (int32_t j=0; j!=num_lhs; j++)
            {
                row[j] = kernel(i,j);
            }

            return row;

        }


        template <class T>
        T* get_kernel_matrix(int32_t &m, int32_t &n, T* target)
        {
            T* result = NULL;

            if (!has_features())
                SG_ERROR( "no features assigned to kernel\n");

            if (target && (m!=get_num_vec_lhs() ||
                        n!=get_num_vec_rhs()) )
            {
                SG_ERROR( "kernel matrix size mismatch\n");
            }

            m=get_num_vec_lhs();
            n=get_num_vec_rhs();

            int64_t total_num = int64_t(m)*n;

            // if lhs == rhs and sizes match assume k(i,j)=k(j,i)
            bool symmetric= (lhs && lhs==rhs && m==n);

            SG_DEBUG( "returning kernel matrix of size %dx%d\n", m, n);

            if (target)
                result=target;
            else
                result=new T[total_num];

            int32_t num_threads=parallel->get_num_threads();
            if (num_threads < 2)
            {
                K_THREAD_PARAM<T> params;
                params.kernel=this;
                params.result=result;
                params.start=0;
                params.end=m;
                params.total_start=0;
                params.total_end=total_num;
                params.n=n;
                params.m=m;
                params.symmetric=symmetric;
                params.verbose=true;
                get_kernel_matrix_helper<T>((void*) &params);
            }
            else
            {
                pthread_t* threads = new pthread_t[num_threads-1];
                K_THREAD_PARAM<T>* params = new K_THREAD_PARAM<T>[num_threads];
                int64_t step= total_num/num_threads;

                int32_t t;

                for (t=0; t<num_threads-1; t++)
                {
                    params[t].kernel = this;
                    params[t].result = result;
                    params[t].start = compute_row_start(t*step, n, symmetric);
                    params[t].end = compute_row_start((t+1)*step, n, symmetric);
                    params[t].total_start=t*step;
                    params[t].total_end=(t+1)*step;
                    params[t].n=n;
                    params[t].m=m;
                    params[t].symmetric=symmetric;
                    params[t].verbose=false;
                    pthread_create(&threads[t], NULL,
                            CKernel::get_kernel_matrix_helper<T>, (void*)&params[t]);
                }

                params[t].kernel = this;
                params[t].result = result;
                params[t].start = compute_row_start(t*step, n, symmetric);
                params[t].end = m;
                params[t].total_start=t*step;
                params[t].total_end=total_num;
                params[t].n=n;
                params[t].m=m;
                params[t].symmetric=symmetric;
                params[t].verbose=true;
                get_kernel_matrix_helper<T>(&params[t]);

                for (t=0; t<num_threads-1; t++)
                    pthread_join(threads[t], NULL);

                delete[] params;
                delete[] threads;
            }

            SG_DONE();

            return result;
        }


        virtual bool init(CFeatures* lhs, CFeatures* rhs);

        virtual bool set_normalizer(CKernelNormalizer* normalizer);

        virtual CKernelNormalizer* get_normalizer();

        virtual bool init_normalizer();

        virtual void cleanup();

        void load(CFile* loader);

        void save(CFile* writer);

        inline CFeatures* get_lhs() { SG_REF(lhs); return lhs; }

        inline CFeatures* get_rhs() { SG_REF(rhs); return rhs; }

        virtual inline int32_t get_num_vec_lhs()
        {
            return num_lhs;
        }

        virtual inline int32_t get_num_vec_rhs()
        {
            return num_rhs;
        }

        virtual inline bool has_features()
        {
            return lhs && rhs;
        }

        inline bool get_lhs_equals_rhs()
        {
            return lhs_equals_rhs;
        }

        virtual void remove_lhs_and_rhs();

        virtual void remove_lhs();

        virtual void remove_rhs();

        virtual EKernelType get_kernel_type()=0 ;

        virtual EFeatureType get_feature_type()=0;

        virtual EFeatureClass get_feature_class()=0;

        inline void set_cache_size(int32_t size)
        {
            cache_size = size;
#ifdef USE_SVMLIGHT
            cache_reset();
#endif //USE_SVMLIGHT
        }

        inline int32_t get_cache_size() { return cache_size; }

#ifdef USE_SVMLIGHT

        inline void cache_reset() { resize_kernel_cache(cache_size); }

        inline int32_t get_max_elems_cache() { return kernel_cache.max_elems; }

        inline int32_t get_activenum_cache() { return kernel_cache.activenum; }

        void get_kernel_row(
            int32_t docnum, int32_t *active2dnum, float64_t *buffer,
            bool full_line=false);

        void cache_kernel_row(int32_t x);

        void cache_multiple_kernel_rows(int32_t* key, int32_t varnum);

        void kernel_cache_reset_lru();

        void kernel_cache_shrink(
            int32_t totdoc, int32_t num_shrink, int32_t *after);

        void resize_kernel_cache(KERNELCACHE_IDX size,
            bool regression_hack=false);

        inline void set_time(int32_t t)
        {
            kernel_cache.time=t;
        }

        inline int32_t kernel_cache_touch(int32_t cacheidx)
        {
            if(kernel_cache.index[cacheidx] != -1)
            {
                kernel_cache.lru[kernel_cache.index[cacheidx]]=kernel_cache.time;
                return(1);
            }
            return(0);
        }

        inline int32_t kernel_cache_check(int32_t cacheidx)
        {
            return(kernel_cache.index[cacheidx] >= 0);
        }

        inline int32_t kernel_cache_space_available()
        {
            return(kernel_cache.elems < kernel_cache.max_elems);
        }

        void kernel_cache_init(int32_t size, bool regression_hack=false);

        void kernel_cache_cleanup();

#endif //USE_SVMLIGHT

        void list_kernel();

        inline bool has_property(EKernelProperty p) { return (properties & p) != 0; }

        virtual void clear_normal();

        virtual void add_to_normal(int32_t vector_idx, float64_t weight);

        inline EOptimizationType get_optimization_type() { return opt_type; }

        virtual inline void set_optimization_type(EOptimizationType t) { opt_type=t;}

        inline bool get_is_initialized() { return optimization_initialized; }

        virtual bool init_optimization(
            int32_t count, int32_t *IDX, float64_t *weights);

        virtual bool delete_optimization();

        bool init_optimization_svm(CSVM * svm) ;

        virtual float64_t compute_optimized(int32_t vector_idx);

        virtual void compute_batch(
            int32_t num_vec, int32_t* vec_idx, float64_t* target,
            int32_t num_suppvec, int32_t* IDX, float64_t* alphas,
            float64_t factor=1.0);

        inline float64_t get_combined_kernel_weight() { return combined_kernel_weight; }

        inline void set_combined_kernel_weight(float64_t nw) { combined_kernel_weight=nw; }

        virtual int32_t get_num_subkernels();

        virtual void compute_by_subkernel(
            int32_t vector_idx, float64_t * subkernel_contrib);

        virtual const float64_t* get_subkernel_weights(int32_t& num_weights);

        virtual void set_subkernel_weights(
            float64_t* weights, int32_t num_weights);

    protected:
        inline void set_property(EKernelProperty p)
        {
            properties |= p;
        }

        inline void unset_property(EKernelProperty p)
        {
            properties &= (properties | p) ^ p;
        }

        inline void set_is_initialized(bool p_init) { optimization_initialized=p_init; }

        virtual float64_t compute(int32_t x, int32_t y)=0;

        int32_t compute_row_start(int64_t offs, int32_t n, bool symmetric)
        {
            int32_t i_start;

            if (symmetric)
                i_start=(int32_t) CMath::floor(n-CMath::sqrt(CMath::sq((float64_t) n)-offs));
            else
                i_start=(int32_t) (offs/int64_t(n));

            return i_start;
        }

        template <class T>
        static void* get_kernel_matrix_helper(void* p)
        {
            K_THREAD_PARAM<T>* params= (K_THREAD_PARAM<T>*) p;
            int32_t i_start=params->start;
            int32_t i_end=params->end;
            CKernel* k=params->kernel;
            T* result=params->result;
            bool symmetric=params->symmetric;
            int32_t n=params->n;
            int32_t m=params->m;
            bool verbose=params->verbose;
            int64_t total_start=params->total_start;
            int64_t total_end=params->total_end;
            int64_t total=total_start;

            for (int32_t i=i_start; i<i_end; i++)
            {
                int32_t j_start=0;

                if (symmetric)
                    j_start=i;

                for (int32_t j=j_start; j<n; j++)
                {
                    float64_t v=k->kernel(i,j);
                    result[i+j*m]=v;

                    if (symmetric && i!=j)
                        result[j+i*m]=v;

                    if (verbose)
                    {
                        total++;

                        if (symmetric && i!=j)
                            total++;

                        if (total%100 == 0)
                            k->SG_PROGRESS(total, total_start, total_end);

                        if (CSignal::cancel_computations())
                            break;
                    }
                }

            }

            return NULL;
        }

        virtual void load_serializable_post() throw (ShogunException);

        virtual void save_serializable_pre() throw (ShogunException);

        virtual void save_serializable_post() throw (ShogunException);

    private:
        void init();


#ifdef USE_SVMLIGHT

#ifndef DOXYGEN_SHOULD_SKIP_THIS
        struct KERNEL_CACHE {
            int32_t   *index;
            int32_t   *invindex;
            int32_t   *active2totdoc;
            int32_t   *totdoc2active;
            int32_t   *lru;
            int32_t   *occu;
            int32_t   elems;
            int32_t   max_elems;
            int32_t   time;
            int32_t   activenum;

            KERNELCACHE_ELEM  *buffer;
            KERNELCACHE_IDX   buffsize;
        };

        struct S_KTHREAD_PARAM
        {
            CKernel* kernel;
            KERNEL_CACHE* kernel_cache;
            KERNELCACHE_ELEM** cache;
            int32_t* uncached_rows;
            int32_t num_uncached;
            uint8_t* needs_computation;
            int32_t start;
            int32_t end;
            int32_t num_vectors;
        };
#endif // DOXYGEN_SHOULD_SKIP_THIS

        static void* cache_multiple_kernel_row_helper(void* p);

        void   kernel_cache_free(int32_t cacheidx);
        int32_t   kernel_cache_malloc();
        int32_t   kernel_cache_free_lru();
        KERNELCACHE_ELEM *kernel_cache_clean_and_malloc(int32_t cacheidx);
#endif //USE_SVMLIGHT


    protected:
        int32_t cache_size;

#ifdef USE_SVMLIGHT

        KERNEL_CACHE kernel_cache;
#endif //USE_SVMLIGHT

        KERNELCACHE_ELEM* kernel_matrix;

        CFeatures* lhs;
        CFeatures* rhs;

        bool lhs_equals_rhs;

        int32_t num_lhs;
        int32_t num_rhs;

        float64_t combined_kernel_weight;

        bool optimization_initialized;
        EOptimizationType opt_type;

        uint64_t  properties;

        CKernelNormalizer* normalizer;
};

}
#endif /* _KERNEL_H__ */