SHOGUN v0.9.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) 2009 Sebastian J. Schultheiss and Soeren Sonnenburg 00008 * Copyright (C) 2009 Max-Planck-Society 00009 */ 00010 00011 #include "lib/common.h" 00012 #include "kernel/RegulatoryModulesStringKernel.h" 00013 #include "features/Features.h" 00014 #include "features/SimpleFeatures.h" 00015 #include "lib/io.h" 00016 00017 using namespace shogun; 00018 00019 CRegulatoryModulesStringKernel::CRegulatoryModulesStringKernel(void) 00020 : CStringKernel<char>(0), width(0.0), degree(0), shift(0), window(0), 00021 motif_positions_lhs(NULL), motif_positions_rhs(NULL), 00022 position_weights(NULL), weights(NULL) 00023 { 00024 SG_UNSTABLE("CRegulatoryModulesStringKernel::" 00025 "CRegulatoryModulesStringKernel(void)", "\n"); 00026 } 00027 00028 CRegulatoryModulesStringKernel::CRegulatoryModulesStringKernel( 00029 int32_t size, float64_t w, int32_t d, int32_t s, int32_t wl) 00030 : CStringKernel<char>(size), width(w), degree(d), shift(s), window(wl), 00031 motif_positions_lhs(NULL), motif_positions_rhs(NULL), position_weights(NULL), weights(NULL) 00032 { 00033 } 00034 00035 CRegulatoryModulesStringKernel::CRegulatoryModulesStringKernel(CStringFeatures<char>* lstr, CStringFeatures<char>* rstr, 00036 CSimpleFeatures<uint16_t>* lpos, CSimpleFeatures<uint16_t>* rpos, 00037 float64_t w, int32_t d, int32_t s, int32_t wl, int32_t size) 00038 : CStringKernel<char>(size), width(w), degree(d), shift(s), window(wl), 00039 motif_positions_lhs(NULL), motif_positions_rhs(NULL), position_weights(NULL), weights(NULL) 00040 { 00041 set_motif_positions(lpos, rpos); 00042 init(lstr,rstr); 00043 } 00044 00045 CRegulatoryModulesStringKernel::~CRegulatoryModulesStringKernel() 00046 { 00047 SG_UNREF(motif_positions_lhs); 00048 SG_UNREF(motif_positions_rhs); 00049 } 00050 00051 bool CRegulatoryModulesStringKernel::init(CFeatures* l, CFeatures* r) 00052 { 00053 ASSERT(motif_positions_lhs); 00054 ASSERT(motif_positions_rhs); 00055 00056 if (l->get_num_vectors() != motif_positions_lhs->get_num_vectors()) 00057 SG_ERROR("Number of vectors does not agree (LHS: %d, Motif LHS: %d).\n", 00058 l->get_num_vectors(), motif_positions_lhs->get_num_vectors()); 00059 if (r->get_num_vectors() != motif_positions_rhs->get_num_vectors()) 00060 SG_ERROR("Number of vectors does not agree (RHS: %d, Motif RHS: %d).\n", 00061 r->get_num_vectors(), motif_positions_rhs->get_num_vectors()); 00062 00063 set_wd_weights(); 00064 CStringKernel<char>::init(l, r); 00065 return init_normalizer(); 00066 } 00067 00068 void CRegulatoryModulesStringKernel::set_motif_positions( 00069 CSimpleFeatures<uint16_t>* positions_lhs, CSimpleFeatures<uint16_t>* positions_rhs) 00070 { 00071 ASSERT(positions_lhs); 00072 ASSERT(positions_rhs); 00073 SG_UNREF(motif_positions_lhs); 00074 SG_UNREF(motif_positions_rhs); 00075 if (positions_lhs->get_num_features() != positions_rhs->get_num_features()) 00076 SG_ERROR("Number of dimensions does not agree.\n"); 00077 00078 motif_positions_lhs=positions_lhs; 00079 motif_positions_rhs=positions_rhs; 00080 SG_REF(positions_lhs); 00081 SG_REF(positions_rhs); 00082 } 00083 00084 float64_t CRegulatoryModulesStringKernel::compute(int32_t idx_a, int32_t idx_b) 00085 { 00086 ASSERT(motif_positions_lhs); 00087 ASSERT(motif_positions_rhs); 00088 00089 int32_t alen, blen; 00090 bool free_avec, free_bvec; 00091 char* avec=((CStringFeatures<char>*) lhs)->get_feature_vector(idx_a, alen, free_avec); 00092 char* bvec=((CStringFeatures<char>*) rhs)->get_feature_vector(idx_b, blen, free_bvec); 00093 00094 int32_t alen_pos, blen_pos; 00095 bool afree_pos, bfree_pos; 00096 uint16_t* positions_a = motif_positions_lhs->get_feature_vector(idx_a, alen_pos, afree_pos); 00097 uint16_t* positions_b = motif_positions_rhs->get_feature_vector(idx_b, blen_pos, bfree_pos); 00098 ASSERT(alen_pos==blen_pos); 00099 int32_t num_pos=alen_pos; 00100 00101 00102 float64_t result_rbf=0; 00103 float64_t result_wds=0; 00104 00105 for (int32_t p=0; p<num_pos; p++) 00106 { 00107 result_rbf+=CMath::sq(positions_a[p]-positions_b[p]); 00108 00109 for (int32_t p2=0; p2<num_pos; p2++) //p+1 and below * 2 00110 result_rbf+=CMath::sq( (positions_a[p]-positions_a[p2]) - (positions_b[p]-positions_b[p2]) ); 00111 00112 int32_t limit = window; 00113 if (window + positions_a[p] > alen) 00114 limit = alen - positions_a[p]; 00115 00116 if (window + positions_b[p] > blen) 00117 limit = CMath::min(limit, blen - positions_b[p]); 00118 00119 result_wds+=compute_wds(&avec[positions_a[p]], &bvec[positions_b[p]], 00120 limit); 00121 } 00122 00123 float64_t result=exp(-result_rbf/width)+result_wds; 00124 00125 ((CStringFeatures<char>*) lhs)->free_feature_vector(avec, idx_a, free_avec); 00126 ((CStringFeatures<char>*) rhs)->free_feature_vector(bvec, idx_b, free_bvec); 00127 ((CSimpleFeatures<uint16_t>*) lhs)->free_feature_vector(positions_a, idx_a, afree_pos); 00128 ((CSimpleFeatures<uint16_t>*) rhs)->free_feature_vector(positions_b, idx_b, bfree_pos); 00129 00130 return result; 00131 } 00132 00133 float64_t CRegulatoryModulesStringKernel::compute_wds( 00134 char* avec, char* bvec, int32_t len) 00135 { 00136 float64_t* max_shift_vec = new float64_t[shift]; 00137 float64_t sum0=0 ; 00138 for (int32_t i=0; i<shift; i++) 00139 max_shift_vec[i]=0 ; 00140 00141 // no shift 00142 for (int32_t i=0; i<len; i++) 00143 { 00144 if ((position_weights!=NULL) && (position_weights[i]==0.0)) 00145 continue ; 00146 00147 float64_t sumi = 0.0 ; 00148 for (int32_t j=0; (j<degree) && (i+j<len); j++) 00149 { 00150 if (avec[i+j]!=bvec[i+j]) 00151 break ; 00152 sumi += weights[j]; 00153 } 00154 if (position_weights!=NULL) 00155 sum0 += position_weights[i]*sumi ; 00156 else 00157 sum0 += sumi ; 00158 } ; 00159 00160 for (int32_t i=0; i<len; i++) 00161 { 00162 for (int32_t k=1; (k<=shift) && (i+k<len); k++) 00163 { 00164 if ((position_weights!=NULL) && (position_weights[i]==0.0) && (position_weights[i+k]==0.0)) 00165 continue ; 00166 00167 float64_t sumi1 = 0.0 ; 00168 // shift in sequence a 00169 for (int32_t j=0; (j<degree) && (i+j+k<len); j++) 00170 { 00171 if (avec[i+j+k]!=bvec[i+j]) 00172 break ; 00173 sumi1 += weights[j]; 00174 } 00175 float64_t sumi2 = 0.0 ; 00176 // shift in sequence b 00177 for (int32_t j=0; (j<degree) && (i+j+k<len); j++) 00178 { 00179 if (avec[i+j]!=bvec[i+j+k]) 00180 break ; 00181 sumi2 += weights[j]; 00182 } 00183 if (position_weights!=NULL) 00184 max_shift_vec[k-1] += position_weights[i]*sumi1 + position_weights[i+k]*sumi2 ; 00185 else 00186 max_shift_vec[k-1] += sumi1 + sumi2 ; 00187 } ; 00188 } 00189 00190 float64_t result = sum0 ; 00191 for (int32_t i=0; i<shift; i++) 00192 result += max_shift_vec[i]/(2*(i+1)) ; 00193 00194 delete[] max_shift_vec; 00195 return result ; 00196 } 00197 00198 void CRegulatoryModulesStringKernel::set_wd_weights() 00199 { 00200 ASSERT(degree>0); 00201 00202 delete[] weights; 00203 weights=new float64_t[degree]; 00204 00205 int32_t i; 00206 float64_t sum=0; 00207 for (i=0; i<degree; i++) 00208 { 00209 weights[i]=degree-i; 00210 sum+=weights[i]; 00211 } 00212 00213 for (i=0; i<degree; i++) 00214 weights[i]/=sum; 00215 }