SHOGUN
v3.2.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 register_params(); 00020 } 00021 00022 CANOVAKernel::CANOVAKernel(int32_t cache, int32_t d) 00023 : CDotKernel(cache), cardinality(d) 00024 { 00025 register_params(); 00026 } 00027 00028 CANOVAKernel::CANOVAKernel( 00029 CDenseFeatures<float64_t>* l, CDenseFeatures<float64_t>* r, int32_t d, int32_t cache) 00030 : CDotKernel(cache), cardinality(d) 00031 { 00032 register_params(); 00033 init(l, r); 00034 } 00035 00036 CANOVAKernel::~CANOVAKernel() 00037 { 00038 cleanup(); 00039 } 00040 00041 bool CANOVAKernel::init(CFeatures* l, CFeatures* r) 00042 { 00043 cleanup(); 00044 00045 bool result = CDotKernel::init(l,r); 00046 00047 init_normalizer(); 00048 return result; 00049 } 00050 00051 float64_t CANOVAKernel::compute(int32_t idx_a, int32_t idx_b) 00052 { 00053 return compute_rec1(idx_a, idx_b); 00054 } 00055 00056 float64_t CANOVAKernel::compute_rec1(int32_t idx_a, int32_t idx_b) 00057 { 00058 int32_t alen, blen; 00059 bool afree, bfree; 00060 00061 float64_t* avec= 00062 ((CDenseFeatures<float64_t>*) lhs)->get_feature_vector(idx_a, alen, afree); 00063 float64_t* bvec= 00064 ((CDenseFeatures<float64_t>*) rhs)->get_feature_vector(idx_b, blen, bfree); 00065 ASSERT(alen==blen) 00066 00067 float64_t result = compute_recursive1(avec, bvec, alen); 00068 00069 ((CDenseFeatures<float64_t>*) lhs)->free_feature_vector(avec, idx_a, afree); 00070 ((CDenseFeatures<float64_t>*) rhs)->free_feature_vector(bvec, idx_b, bfree); 00071 00072 return result; 00073 } 00074 00075 float64_t CANOVAKernel::compute_rec2(int32_t idx_a, int32_t idx_b) 00076 { 00077 int32_t alen, blen; 00078 bool afree, bfree; 00079 00080 float64_t* avec= 00081 ((CDenseFeatures<float64_t>*) lhs)->get_feature_vector(idx_a, alen, afree); 00082 float64_t* bvec= 00083 ((CDenseFeatures<float64_t>*) rhs)->get_feature_vector(idx_b, blen, bfree); 00084 ASSERT(alen==blen) 00085 00086 float64_t result = compute_recursive2(avec, bvec, alen); 00087 00088 ((CDenseFeatures<float64_t>*) lhs)->free_feature_vector(avec, idx_a, afree); 00089 ((CDenseFeatures<float64_t>*) rhs)->free_feature_vector(bvec, idx_b, bfree); 00090 00091 return result; 00092 } 00093 00094 void CANOVAKernel::register_params() 00095 { 00096 SG_ADD(&cardinality, "cardinality", "Kernel cardinality.", MS_AVAILABLE); 00097 } 00098 00099 00100 float64_t CANOVAKernel::compute_recursive1(float64_t* avec, float64_t* bvec, int32_t len) 00101 { 00102 int32_t DP_len=(cardinality+1)*(len+1); 00103 float64_t* DP = SG_MALLOC(float64_t, DP_len); 00104 00105 ASSERT(DP) 00106 int32_t d=cardinality; 00107 int32_t offs=cardinality+1; 00108 00109 ASSERT(DP_len==(len+1)*offs) 00110 00111 for (int32_t j=0; j < len+1; j++) 00112 DP[j] = 1.0; 00113 00114 for (int32_t k=1; k < d+1; k++) 00115 { 00116 // TRAP d>len case 00117 if (k-1>=len) 00118 return 0.0; 00119 00120 DP[k*offs+k-1] = 0; 00121 for (int32_t j=k; j < len+1; j++) 00122 DP[k*offs+j]=DP[k*offs+j-1]+avec[j-1]*bvec[j-1]*DP[(k-1)*offs+j-1]; 00123 } 00124 00125 float64_t result=DP[d*offs+len]; 00126 00127 SG_FREE(DP); 00128 00129 return result; 00130 } 00131 00132 float64_t CANOVAKernel::compute_recursive2(float64_t* avec, float64_t* bvec, int32_t len) 00133 { 00134 float64_t* KD = SG_MALLOC(float64_t, cardinality+1); 00135 float64_t* KS = SG_MALLOC(float64_t, cardinality+1); 00136 float64_t* vec_pow = SG_MALLOC(float64_t, len); 00137 00138 ASSERT(vec_pow) 00139 ASSERT(KS) 00140 ASSERT(KD) 00141 00142 int32_t d=cardinality; 00143 for (int32_t i=0; i < len; i++) 00144 vec_pow[i] = 1; 00145 00146 for (int32_t k=1; k < d+1; k++) 00147 { 00148 KS[k] = 0; 00149 for (int32_t i=0; i < len; i++) 00150 { 00151 vec_pow[i] *= avec[i]*bvec[i]; 00152 KS[k] += vec_pow[i]; 00153 } 00154 } 00155 00156 KD[0] = 1; 00157 for (int32_t k=1; k < d+1; k++) 00158 { 00159 float64_t sum = 0; 00160 for (int32_t s=1; s < k+1; s++) 00161 { 00162 float64_t sign = 1.0; 00163 if (s % 2 == 0) 00164 sign = -1.0; 00165 00166 sum += sign*KD[k-s]*KS[s]; 00167 } 00168 00169 KD[k] = sum / k; 00170 } 00171 float64_t result=KD[d]; 00172 SG_FREE(vec_pow); 00173 SG_FREE(KS); 00174 SG_FREE(KD); 00175 00176 return result; 00177 }