00001
00002
00003
00004
00005
00006
00007
00008
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
00118 DP_len=(cardinality+1)*(num_feat+1);
00119 DP = SG_MALLOC(float64_t, DP_len);
00120
00121
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
00130 SG_FREE(DP);
00131 DP=NULL;
00132 DP_len=0;
00133
00134
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
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 }