diff --git a/src/interfaces/modular/Kernel.i b/src/interfaces/modular/Kernel.i index 43d80536d64..10879ee1632 100644 --- a/src/interfaces/modular/Kernel.i +++ b/src/interfaces/modular/Kernel.i @@ -59,6 +59,7 @@ PROTOCOLS_CUSTOMKERNEL(CustomKernel, float32_t, "f\0", NPY_FLOAT32) %rename(DistanceKernel) CDistanceKernel; %rename(FixedDegreeStringKernel) CFixedDegreeStringKernel; %rename(GaussianKernel) CGaussianKernel; +%rename(ShiftInvariantKernel) CShiftInvariantKernel; %rename(GaussianCompactKernel) CGaussianCompactKernel; %rename(DirectorKernel) CDirectorKernel; %rename(WaveletKernel) CWaveletKernel; @@ -187,6 +188,7 @@ namespace shogun %include %include %include +%include %include %include %include diff --git a/src/interfaces/modular/Kernel_includes.i b/src/interfaces/modular/Kernel_includes.i index 47efb8f04f0..c532914e8e3 100644 --- a/src/interfaces/modular/Kernel_includes.i +++ b/src/interfaces/modular/Kernel_includes.i @@ -26,6 +26,7 @@ #include #include #include +#include #include #include #include diff --git a/src/shogun/kernel/GaussianCompactKernel.cpp b/src/shogun/kernel/GaussianCompactKernel.cpp index 53f7abf5b07..8e813ce3922 100644 --- a/src/shogun/kernel/GaussianCompactKernel.cpp +++ b/src/shogun/kernel/GaussianCompactKernel.cpp @@ -1,4 +1,5 @@ #include +#include using namespace shogun; diff --git a/src/shogun/kernel/GaussianKernel.cpp b/src/shogun/kernel/GaussianKernel.cpp index 1530d66fae1..80dd2c636a3 100644 --- a/src/shogun/kernel/GaussianKernel.cpp +++ b/src/shogun/kernel/GaussianKernel.cpp @@ -7,6 +7,7 @@ * Written (W) 1999-2010 Soeren Sonnenburg * Written (W) 2011 Abhinav Maurya * Written (W) 2012 Heiko Strathmann + * Written (W) 2016 Soumyajit De * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society * Copyright (C) 2010 Berlin Institute of Technology */ @@ -15,45 +16,36 @@ #include #include #include +#include #include #include using namespace shogun; -void CGaussianKernel::set_width(float64_t w) +CGaussianKernel::CGaussianKernel() : CShiftInvariantKernel() { - REQUIRE(w>0, "width (%f) must be positive\n",w); - m_log_width=CMath::log(w/2.0)/2.0; + register_params(); } -float64_t CGaussianKernel::get_width() const +CGaussianKernel::CGaussianKernel(float64_t w) : CShiftInvariantKernel() { - return CMath::exp(m_log_width*2.0)*2.0; -} - -CGaussianKernel::CGaussianKernel() : CDotKernel() -{ - init(); -} - -CGaussianKernel::CGaussianKernel(float64_t w) : CDotKernel() -{ - init(); + register_params(); set_width(w); } -CGaussianKernel::CGaussianKernel(int32_t size, float64_t w) : CDotKernel(size) +CGaussianKernel::CGaussianKernel(int32_t size, float64_t w) : CShiftInvariantKernel() { - init(); + register_params(); + set_cache_size(size); set_width(w); } -CGaussianKernel::CGaussianKernel(CDotFeatures* l, CDotFeatures* r, - float64_t w, int32_t size) : CDotKernel(size) +CGaussianKernel::CGaussianKernel(CDotFeatures* l, CDotFeatures* r, float64_t w, int32_t size) : CShiftInvariantKernel(l, r) { - init(); + register_params(); + set_cache_size(size); set_width(w); - init(l,r); + init(l, r); } CGaussianKernel::~CGaussianKernel() @@ -74,94 +66,68 @@ CGaussianKernel* CGaussianKernel::obtain_from_generic(CKernel* kernel) return (CGaussianKernel*)kernel; } -#include -CSGObject *CGaussianKernel::shallow_copy() const -{ - // TODO: remove this after all the classes get shallow_copy properly implemented - // this assert is to avoid any subclass of CGaussianKernel accidentally called - // with the implement here - ASSERT(typeid(*this) == typeid(CGaussianKernel)) - CGaussianKernel *ker = new CGaussianKernel(cache_size, get_width()); - if (lhs) - { - ker->init(lhs, rhs); - } - return ker; -} - void CGaussianKernel::cleanup() { - if (sq_lhs != sq_rhs) - SG_FREE(sq_rhs); - sq_rhs = NULL; - - SG_FREE(sq_lhs); - sq_lhs = NULL; - + m_distance->reset_precompute(); CKernel::cleanup(); } -void CGaussianKernel::precompute_squared_helper(float64_t* &buf, CDotFeatures* df) +bool CGaussianKernel::init(CFeatures* l, CFeatures* r) { - ASSERT(df) - int32_t num_vec=df->get_num_vectors(); - buf=SG_MALLOC(float64_t, num_vec); + CShiftInvariantKernel::init(l, r); + m_distance->reset_precompute(); - for (int32_t i=0; idot(i,df, i); -} + REQUIRE(l->has_property(FP_DOT), "Left hand side (%s) must be a subclass of DotFeatures!\n", l->get_name()); + REQUIRE(r->has_property(FP_DOT), "Right hand side (%s) must be a subclass of DotFeatures!\n", r->get_name()) -bool CGaussianKernel::init(CFeatures* l, CFeatures* r) -{ - ///free sq_{r,l}hs first - cleanup(); + int32_t lhs_dim_feature_space=static_cast(l)->get_dim_feature_space(); + int32_t rhs_dim_feature_space=static_cast(r)->get_dim_feature_space(); - CDotKernel::init(l, r); - precompute_squared(); + REQUIRE(lhs_dim_feature_space==rhs_dim_feature_space, + "Train or test features #dimension mismatch (l:%d vs. r:%d)\n", + lhs_dim_feature_space, rhs_dim_feature_space); + + precompute_squared_norms(); return init_normalizer(); } float64_t CGaussianKernel::compute(int32_t idx_a, int32_t idx_b) { - float64_t result=distance(idx_a,idx_b); + float64_t result=distance(idx_a, idx_b); return CMath::exp(-result); } void CGaussianKernel::load_serializable_post() throw (ShogunException) { CKernel::load_serializable_post(); - precompute_squared(); + precompute_squared_norms(); } -void CGaussianKernel::precompute_squared() +void CGaussianKernel::precompute_squared_norms() { - if (!lhs || !rhs) - return; - - precompute_squared_helper(sq_lhs, (CDotFeatures*) lhs); - - if (lhs==rhs) - sq_rhs=sq_lhs; - else - precompute_squared_helper(sq_rhs, (CDotFeatures*) rhs); + if (lhs && rhs) + { + m_distance->precompute_lhs(); + m_distance->precompute_rhs(); + } } -SGMatrix CGaussianKernel::get_parameter_gradient( - const TParameter* param, index_t index) +SGMatrix CGaussianKernel::get_parameter_gradient(const TParameter* param, index_t index) { - REQUIRE(lhs && rhs, "Features not set!\n") + REQUIRE(lhs, "The left hand side feature instance cannot be NULL!\n"); + REQUIRE(rhs, "The right hand side feature instance cannot be NULL!\n"); if (!strcmp(param->m_name, "log_width")) { SGMatrix derivative=SGMatrix(num_lhs, num_rhs); - for (int j=0; j CGaussianKernel::get_parameter_gradient( } } -void CGaussianKernel::init() +void CGaussianKernel::register_params() { set_width(1.0); - sq_lhs=NULL; - sq_rhs=NULL; + set_cache_size(10); + m_distance=new CEuclideanDistance(); + SG_REF(m_distance); SG_ADD(&m_log_width, "log_width", "Kernel width in log domain", MS_AVAILABLE, GRADIENT_AVAILABLE); } -float64_t CGaussianKernel::distance(int32_t idx_a, int32_t idx_b) +void CGaussianKernel::set_width(float64_t w) { - return (sq_lhs[idx_a]+sq_rhs[idx_b]-2*CDotKernel::compute(idx_a,idx_b))/get_width(); + REQUIRE(w>0, "width (%f) must be positive\n",w); + m_log_width=CMath::log(w/2.0)/2.0; +} + +float64_t CGaussianKernel::get_width() const +{ + return CMath::exp(m_log_width*2.0)*2.0; +} + +float64_t CGaussianKernel::distance(int32_t idx_a, int32_t idx_b) const +{ + float64_t distance=CShiftInvariantKernel::distance(idx_a, idx_b); + return distance*distance/get_width(); +} + +#include +CSGObject *CGaussianKernel::shallow_copy() const +{ + // TODO: remove this after all the classes get shallow_copy properly implemented + // this assert is to avoid any subclass of CGaussianKernel accidentally called + // with the implement here + ASSERT(typeid(*this) == typeid(CGaussianKernel)) + CGaussianKernel *ker = new CGaussianKernel(cache_size, get_width()); + if (lhs) + ker->init(lhs, rhs); + + return ker; } diff --git a/src/shogun/kernel/GaussianKernel.h b/src/shogun/kernel/GaussianKernel.h index 4c9c21f0c0f..3a177cc697a 100644 --- a/src/shogun/kernel/GaussianKernel.h +++ b/src/shogun/kernel/GaussianKernel.h @@ -7,6 +7,7 @@ * Written (W) 1999-2010 Soeren Sonnenburg * Written (W) 2011 Abhinav Maurya * Written (W) 2012 Heiko Strathmann + * Written (W) 2016 Soumyajit De * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society * Copyright (C) 2010 Berlin Institute of Technology */ @@ -15,15 +16,16 @@ #define GAUSSIANKERNEL_H #include - #include -#include -#include -#include +#include +#include namespace shogun { - class CDotFeatures; + +class CFeatures; +class CDotFeatures; + /** @brief The well known Gaussian kernel (swiss army knife for SVMs) computed * on CDotFeatures. * @@ -36,149 +38,166 @@ namespace shogun * where \f$\tau\f$ is the kernel width. * */ -class CGaussianKernel: public CDotKernel +class CGaussianKernel: public CShiftInvariantKernel { - public: - /** default constructor */ - CGaussianKernel(); - - /** constructor - * - * @param width width - */ - CGaussianKernel(float64_t width); - - /** constructor - * - * @param size cache size - * @param width width - */ - CGaussianKernel(int32_t size, float64_t width); - - /** constructor - * - * @param l features of left-hand side - * @param r features of right-hand side - * @param width width - * @param size cache size - */ - CGaussianKernel(CDotFeatures* l, CDotFeatures* r, - float64_t width, int32_t size=10); - - virtual ~CGaussianKernel(); - - /** @param kernel is casted to CGaussianKernel, error if not possible - * is SG_REF'ed - * @return casted CGaussianKernel object - */ - static CGaussianKernel* obtain_from_generic(CKernel* kernel); - - /** Make a shallow copy of the kernel */ - virtual CSGObject* shallow_copy() const; - - /** initialize kernel - * - * @param l features of left-hand side - * @param r features of right-hand side - * @return if initializing was successful - */ - virtual bool init(CFeatures* l, CFeatures* r); - - /** clean up kernel */ - virtual void cleanup(); - - /** return what type of kernel we are - * - * @return kernel type GAUSSIAN - */ - virtual EKernelType get_kernel_type() { return K_GAUSSIAN; } - - /** return the kernel's name - * - * @return name Gaussian - */ - virtual const char* get_name() const { return "GaussianKernel"; } - - /** set the kernel's width - * - * @param w kernel width - */ - virtual void set_width(float64_t w); - - /** return the kernel's width - * - * @return kernel width - */ - virtual float64_t get_width() const; - - /** return derivative with respect to specified parameter - * - * @param param the parameter - * @param index the index of the element if parameter is a vector - * - * @return gradient with respect to parameter - */ - virtual SGMatrix get_parameter_gradient( - const TParameter* param, index_t index=-1); - - protected: - /** compute kernel function for features a and b - * idx_{a,b} denote the index of the feature vectors - * in the corresponding feature object - * - * @param idx_a index a - * @param idx_b index b - * @return computed kernel function at indices a,b - */ - virtual float64_t compute(int32_t idx_a, int32_t idx_b); - - /** Can (optionally) be overridden to post-initialize some member - * variables which are not PARAMETER::ADD'ed. Make sure that at first - * the overridden method BASE_CLASS::LOAD_SERIALIZABLE_POST is called. - * - * @exception ShogunException Will be thrown if an error occurres. - */ - virtual void load_serializable_post() throw (ShogunException); - - - /** compute the distance between features a and b - * idx_{a,b} denote the index of the feature vectors - * in the corresponding feature object - * - * @param idx_a index a - * @param idx_b index b - * @return computed the distance - * - * Note that in GaussianKernel, - * kernel(idx_a, idx_b)=exp(-distance(idx_a, idx_b)) - * \f[ - * distance({\bf x},{\bf y})= \frac{||{\bf x}-{\bf y}||^2}{\tau} - * \f] - */ - virtual float64_t distance(int32_t idx_a, int32_t idx_b); - private: - /** helper function to compute quadratic terms in - * (a-b)^2 (== a^2+b^2-2ab) - */ - void precompute_squared(); - - /** helper function to compute quadratic terms in - * (a-b)^2 (== a^2+b^2-2ab) - * - * @param buf buffer to store squared terms (will be allocated) - * @param df dot feature object based on which k(i,i) is computed - * */ - void precompute_squared_helper(float64_t* &buf, CDotFeatures* df); - - void init(); - - protected: - /** width */ - float64_t m_log_width; - /** squared left-hand side */ - float64_t* sq_lhs; - /** squared right-hand side */ - float64_t* sq_rhs; +public: + /** default constructor */ + CGaussianKernel(); + + /** constructor + * + * @param width width + */ + CGaussianKernel(float64_t width); + + /** constructor + * + * @param size cache size + * @param width width + */ + CGaussianKernel(int32_t size, float64_t width); + + /** constructor + * + * @param l features of left-hand side + * @param r features of right-hand side + * @param width width + * @param size cache size + */ + CGaussianKernel(CDotFeatures* l, CDotFeatures* r, float64_t width, int32_t size=10); + + /** destructor */ + virtual ~CGaussianKernel(); + + /** @param kernel is casted to CGaussianKernel, error if not possible + * is SG_REF'ed + * @return casted CGaussianKernel object + */ + static CGaussianKernel* obtain_from_generic(CKernel* kernel); + + /** Make a shallow copy of the kernel */ + virtual CSGObject* shallow_copy() const; + + /** initialize kernel + * + * @param l features of left-hand side + * @param r features of right-hand side + * @return if initializing was successful + */ + virtual bool init(CFeatures* l, CFeatures* r); + + /** clean up kernel */ + virtual void cleanup(); + + /** return what type of kernel we are + * + * @return kernel type GAUSSIAN + */ + virtual EKernelType get_kernel_type() + { + return K_GAUSSIAN; + } + + /** return feature class the kernel can deal with + * + * Caussian kernel returns unknown since features can be based on anything + * + * @return feature class ANY + */ + virtual EFeatureClass get_feature_class() + { + return C_ANY; + } + + /** return feature type the kernel can deal with + * + * Caussian kernel returns unknown since features can be based on anything + * + * @return ANY feature type + */ + virtual EFeatureType get_feature_type() + { + return F_ANY; + } + + /** return the kernel's name + * + * @return name Gaussian + */ + virtual const char* get_name() const + { + return "GaussianKernel"; + } + + /** set the kernel's width + * + * @param w kernel width + */ + virtual void set_width(float64_t w); + + /** return the kernel's width + * + * @return kernel width + */ + virtual float64_t get_width() const; + + /** return derivative with respect to specified parameter + * + * @param param the parameter + * @param index the index of the element if parameter is a vector + * + * @return gradient with respect to parameter + */ + virtual SGMatrix get_parameter_gradient(const TParameter* param, index_t index=-1); + +protected: + /** compute kernel function for features a and b + * idx_{a,b} denote the index of the feature vectors + * in the corresponding feature object + * + * @param idx_a index a + * @param idx_b index b + * @return computed kernel function at indices a,b + */ + virtual float64_t compute(int32_t idx_a, int32_t idx_b); + + /** Can (optionally) be overridden to post-initialize some member + * variables which are not PARAMETER::ADD'ed. Make sure that at first + * the overridden method BASE_CLASS::LOAD_SERIALIZABLE_POST is called. + * + * @exception ShogunException Will be thrown if an error occurres. + */ + virtual void load_serializable_post() throw (ShogunException); + + /** compute the distance between features a and b + * idx_{a,b} denote the index of the feature vectors + * in the corresponding feature object + * + * @param idx_a index a + * @param idx_b index b + * @return computed the distance + * + * Note that in GaussianKernel, + * kernel(idx_a, idx_b)=exp(-distance(idx_a, idx_b)) + * \f[ + * distance({\bf x},{\bf y})= \frac{||{\bf x}-{\bf y}||_{l_2}^2}{\tau} + * \f] + */ + virtual float64_t distance(int32_t idx_a, int32_t idx_b) const; + + /** width */ + float64_t m_log_width; + +private: + /** helper function to compute quadratic terms in (a-b)^2 (== a^2+b^2-2ab), + * internally calls precompute_lhs and precomputed_rhs on the Euclidean + * distance object that is inside. + */ + void precompute_squared_norms(); + + /** initializes with default values and registers parameters */ + void register_params(); + }; } #endif /* _GAUSSIANKERNEL_H__ */ diff --git a/src/shogun/kernel/ShiftInvariantKernel.cpp b/src/shogun/kernel/ShiftInvariantKernel.cpp index 955b62dbc15..398c60380f3 100644 --- a/src/shogun/kernel/ShiftInvariantKernel.cpp +++ b/src/shogun/kernel/ShiftInvariantKernel.cpp @@ -53,9 +53,9 @@ CShiftInvariantKernel::~CShiftInvariantKernel() bool CShiftInvariantKernel::init(CFeatures* l, CFeatures* r) { - REQUIRE(m_distance, "The distance instance cannot be NULL!\n"); CKernel::init(l,r); - m_distance->init(l, r); + if (m_distance) + m_distance->init(l, r); return init_normalizer(); } diff --git a/tests/unit/kernel/GaussianKernel_unittest.cc b/tests/unit/kernel/GaussianKernel_unittest.cc new file mode 100644 index 00000000000..8588e3873d8 --- /dev/null +++ b/tests/unit/kernel/GaussianKernel_unittest.cc @@ -0,0 +1,115 @@ +/* + * Copyright (c) The Shogun Machine Learning Toolbox + * Written (w) 2016 Soumyajit De + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + * The views and conclusions contained in the software and documentation are those + * of the authors and should not be interpreted as representing official policies, + * either expressed or implied, of the Shogun Development Team. + */ + +#include + +#ifdef HAVE_CXX11 + +#include +#include +#include +#include +#include +#include +#include + +using namespace shogun; +using std::iota; +using std::for_each; + +TEST(GaussianKernel, compute) +{ + const index_t dim=1; + const index_t N=4; + const index_t M=6; + const float64_t width=2*pow(2, 0.05); + + SGMatrix data_1(dim, N); + SGMatrix data_2(dim, M); + + iota(data_1.data(), data_1.data()+data_1.size(), 1); + iota(data_2.data(), data_2.data()+data_2.size(), data_1.size()+1); + + for_each(data_1.data(), data_1.data()+data_1.size(), [&data_1](float64_t& val) { val=val/data_1.size(); }); + for_each(data_2.data(), data_2.data()+data_2.size(), [&data_2](float64_t& val) { val=val/data_2.size(); }); + + auto feats_1=some >(data_1); + auto feats_2=some >(data_2); + + auto kernel=some(10, width); + kernel->init(feats_1, feats_2); + + auto euclidean_distance=some(); + euclidean_distance->init(feats_1, feats_2); + + for (auto i=0; idistance(i, j); + EXPECT_NEAR(kernel->kernel(i, j), CMath::exp(-distance*distance/width), 1E-6); + } + } +} + +TEST(GaussianKernel, precomputed_distance) +{ + const index_t dim=1; + const index_t N=4; + const index_t M=6; + const float64_t width=0.05; + + SGMatrix data_1(dim, N); + SGMatrix data_2(dim, M); + + iota(data_1.data(), data_1.data()+data_1.size(), 1); + iota(data_2.data(), data_2.data()+data_2.size(), data_1.size()+1); + + for_each(data_1.data(), data_1.data()+data_1.size(), [&data_1](float64_t& val) { val=val/data_1.size(); }); + for_each(data_2.data(), data_2.data()+data_2.size(), [&data_2](float64_t& val) { val=val/data_2.size(); }); + + auto feats_1=some >(data_1); + auto feats_2=some >(data_2); + + auto kernel_1=some(10, width); + auto kernel_2=some(10, width); + + kernel_1->init(feats_1, feats_2); + kernel_2->init(feats_1, feats_2); + + kernel_1->precompute_distance(); + + for (auto i=0; ikernel(i, j), kernel_2->kernel(i, j), 1E-6); + } +} +#endif // HAVE_CXX11