SHOGUN
v1.1.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * Written (W) 2011 Andrew Tereskin 00008 * Copyright (C) 2011 Berlin Institute of Technology and Max-Planck-Society 00009 */ 00010 00011 #include <math.h> 00012 #include <shogun/kernel/ANOVAKernel.h> 00013 #include <shogun/mathematics/Math.h> 00014 00015 using namespace shogun; 00016 00017 CANOVAKernel::CANOVAKernel(): CDotKernel(0), cardinality(1.0) 00018 { 00019 init(); 00020 register_params(); 00021 } 00022 00023 CANOVAKernel::CANOVAKernel(int32_t cache, int32_t d) 00024 : CDotKernel(cache), cardinality(d) 00025 { 00026 init(); 00027 register_params(); 00028 } 00029 00030 CANOVAKernel::CANOVAKernel( 00031 CSimpleFeatures<float64_t>* l, CSimpleFeatures<float64_t>* r, int32_t d, int32_t cache) 00032 : CDotKernel(cache), cardinality(d) 00033 { 00034 init(); 00035 register_params(); 00036 init(l, r); 00037 } 00038 00039 CANOVAKernel::~CANOVAKernel() 00040 { 00041 cleanup(); 00042 } 00043 00044 bool CANOVAKernel::init(CFeatures* l, CFeatures* r) 00045 { 00046 cleanup(); 00047 00048 bool result = CDotKernel::init(l,r); 00049 00050 allocate_arrays(); 00051 init_normalizer(); 00052 return result; 00053 } 00054 00055 float64_t CANOVAKernel::compute(int32_t idx_a, int32_t idx_b) 00056 { 00057 return compute_rec1(idx_a, idx_b); 00058 } 00059 00060 float64_t CANOVAKernel::compute_rec1(int32_t idx_a, int32_t idx_b) 00061 { 00062 int32_t alen, blen; 00063 bool afree, bfree; 00064 00065 float64_t* avec= 00066 ((CSimpleFeatures<float64_t>*) lhs)->get_feature_vector(idx_a, alen, afree); 00067 float64_t* bvec= 00068 ((CSimpleFeatures<float64_t>*) rhs)->get_feature_vector(idx_b, blen, bfree); 00069 ASSERT(alen==blen); 00070 00071 float64_t result = compute_recursive1(avec, bvec, alen); 00072 00073 ((CSimpleFeatures<float64_t>*) lhs)->free_feature_vector(avec, idx_a, afree); 00074 ((CSimpleFeatures<float64_t>*) rhs)->free_feature_vector(bvec, idx_b, bfree); 00075 00076 return result; 00077 } 00078 00079 float64_t CANOVAKernel::compute_rec2(int32_t idx_a, int32_t idx_b) 00080 { 00081 int32_t alen, blen; 00082 bool afree, bfree; 00083 00084 float64_t* avec= 00085 ((CSimpleFeatures<float64_t>*) lhs)->get_feature_vector(idx_a, alen, afree); 00086 float64_t* bvec= 00087 ((CSimpleFeatures<float64_t>*) rhs)->get_feature_vector(idx_b, blen, bfree); 00088 ASSERT(alen==blen); 00089 00090 float64_t result = compute_recursive2(avec, bvec, alen); 00091 00092 ((CSimpleFeatures<float64_t>*) lhs)->free_feature_vector(avec, idx_a, afree); 00093 ((CSimpleFeatures<float64_t>*) rhs)->free_feature_vector(bvec, idx_b, bfree); 00094 00095 return result; 00096 } 00097 00098 void CANOVAKernel::init() 00099 { 00101 DP=NULL; 00102 00104 KD=NULL; 00105 KS=NULL; 00106 vec_pow=NULL; 00107 } 00108 00109 void CANOVAKernel::allocate_arrays() 00110 { 00111 cleanup(); 00112 00113 ASSERT(lhs && rhs); 00114 int32_t num_feat = ((CSimpleFeatures<float64_t>*) lhs)->get_num_features(); 00115 ASSERT(num_feat == ((CSimpleFeatures<float64_t>*) rhs)->get_num_features()); 00116 00117 //compute_recursive1 00118 DP_len=(cardinality+1)*(num_feat+1); 00119 DP = SG_MALLOC(float64_t, DP_len); 00120 00121 //compute_recursive2 00122 KD = SG_MALLOC(float64_t, cardinality+1); 00123 KS = SG_MALLOC(float64_t, cardinality+1); 00124 vec_pow = SG_MALLOC(float64_t, num_feat); 00125 } 00126 00127 void CANOVAKernel::cleanup() 00128 { 00129 //compute_recursive1 00130 SG_FREE(DP); 00131 DP=NULL; 00132 DP_len=0; 00133 00134 //compute_recursive2 00135 SG_FREE(KD); 00136 KD=NULL; 00137 SG_FREE(KS); 00138 KS=NULL; 00139 SG_FREE(vec_pow); 00140 vec_pow=NULL; 00141 } 00142 00143 void CANOVAKernel::load_serializable_post() throw (ShogunException) 00144 { 00145 CKernel::load_serializable_post(); 00146 allocate_arrays(); 00147 } 00148 00149 void CANOVAKernel::register_params() 00150 { 00151 m_parameters->add(&cardinality, "cardinality", "Kernel cardinality."); 00152 } 00153 00154 00155 float64_t CANOVAKernel::compute_recursive1(float64_t* avec, float64_t* bvec, int32_t len) 00156 { 00157 ASSERT(DP); 00158 int32_t d=cardinality; 00159 int32_t offs=cardinality+1; 00160 00161 ASSERT(DP_len==(len+1)*offs); 00162 00163 for (int32_t j=0; j < len+1; j++) 00164 DP[j] = 1.0; 00165 00166 for (int32_t k=1; k < d+1; k++) 00167 { 00168 // TRAP d>len case 00169 if (k-1>=len) 00170 return 0.0; 00171 00172 DP[k*offs+k-1] = 0; 00173 for (int32_t j=k; j < len+1; j++) 00174 DP[k*offs+j]=DP[k*offs+j-1]+avec[j-1]*bvec[j-1]*DP[(k-1)*offs+j-1]; 00175 } 00176 00177 float64_t result=DP[d*offs+len]; 00178 00179 return result; 00180 } 00181 00182 float64_t CANOVAKernel::compute_recursive2(float64_t* avec, float64_t* bvec, int32_t len) 00183 { 00184 ASSERT(vec_pow); 00185 ASSERT(KS); 00186 ASSERT(KD); 00187 00188 int32_t d=cardinality; 00189 for (int32_t i=0; i < len; i++) 00190 vec_pow[i] = 1; 00191 00192 for (int32_t k=1; k < d+1; k++) 00193 { 00194 KS[k] = 0; 00195 for (int32_t i=0; i < len; i++) 00196 { 00197 vec_pow[i] *= avec[i]*bvec[i]; 00198 KS[k] += vec_pow[i]; 00199 } 00200 } 00201 00202 KD[0] = 1; 00203 for (int32_t k=1; k < d+1; k++) 00204 { 00205 float64_t sum = 0; 00206 for (int32_t s=1; s < k+1; s++) 00207 { 00208 float64_t sign = 1.0; 00209 if (s % 2 == 0) 00210 sign = -1.0; 00211 00212 sum += sign*KD[k-s]*KS[s]; 00213 } 00214 00215 KD[k] = sum / k; 00216 } 00217 float64_t result=KD[d]; 00218 00219 return result; 00220 }