Kernel.h

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) 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 <shogun/lib/common.h>
#include <shogun/lib/Signal.h>
#include <shogun/io/File.h>
#include <shogun/mathematics/Math.h>
#include <shogun/features/FeatureTypes.h>
#include <shogun/base/SGObject.h>
#include <shogun/features/Features.h>
#include <shogun/kernel/KernelNormalizer.h>

#include <vector>

namespace shogun
{
    class CFile;
    class CFeatures;
    class CKernelNormalizer;

#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_SPECTRUMRBF = 122,
    K_SPECTRUMMISMATCHRBF = 123,
    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,
    K_WAVELET = 310,
    K_WAVE = 320,
    K_CAUCHY = 330,
    K_TSTUDENT = 340,
    K_RATIONAL_QUADRATIC = 350,
    K_MULTIQUADRIC = 360,
    K_EXPONENTIAL = 370,
    K_SPHERICAL = 380,
    K_SPLINE = 390,
    K_ANOVA = 400,
    K_POWER = 410,
    K_LOG = 420,
    K_CIRCULAR = 430,
    K_INVERSEMULTIQUADRIC = 440,
    K_DISTANTSEGMENTS = 450,
    K_BESSEL = 460,
};

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
};

#ifndef DOXYGEN_SHOULD_SKIP_THIS

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;
};
#endif

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);
        }

        SGMatrix<float64_t> get_kernel_matrix()
        {
            return get_kernel_matrix<float64_t>();
        }

        virtual SGVector<float64_t> get_kernel_col(int32_t j)
        {

            SGVector<float64_t> col = SGVector<float64_t>(num_rhs);

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

            return col;
        }


        virtual SGVector<float64_t> get_kernel_row(int32_t i)
        {
            SGVector<float64_t> row = SGVector<float64_t>(num_lhs);

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

            return row;
        }

        template <class T>
        SGMatrix<T> get_kernel_matrix()
        {
            T* result = NULL;

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

            int32_t m=get_num_vec_lhs();
            int32_t 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);

            result=SG_MALLOC(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 = SG_MALLOC(pthread_t, num_threads-1);
                K_THREAD_PARAM<T>* params = SG_MALLOC(K_THREAD_PARAM<T>, num_threads);
                int64_t step= total_num/num_threads;

                int32_t t;

                num_threads--;
                for (t=0; t<num_threads; 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;

                    int code=pthread_create(&threads[t], NULL,
                            CKernel::get_kernel_matrix_helper<T>, (void*)&params[t]);

                    if (code != 0)
                    {
                        SG_WARNING("Thread creation failed (thread %d of %d) "
                                "with error:'%s'\n",t, num_threads, strerror(code));
                        num_threads=t;
                        break;
                    }
                }

                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; t++)
                {
                    if (pthread_join(threads[t], NULL) != 0)
                        SG_WARNING("pthread_join of thread %d/%d failed\n", t, num_threads);
                }

                SG_FREE(params);
                SG_FREE(threads);
            }

            SG_DONE();

            return SGMatrix<T>(result,m,n,true);
        }


        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;

        }

        inline int32_t get_cache_size() { return cache_size; }



        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(SGVector<float64_t> 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);
        virtual void register_params();

    private:
        void init();




    protected:
        int32_t cache_size;



        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__ */