SHOGUN  6.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules
SpectrumRBFKernel.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 1999-2009 Soeren Sonnenburg
8  * Written (W) 1999-2008 Gunnar Raetsch
9  * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society
10  */
11 
12 #include <vector>
13 
14 #include <shogun/lib/common.h>
15 #include <shogun/io/SGIO.h>
16 #include <shogun/lib/Signal.h>
17 #include <shogun/lib/Trie.h>
18 #include <shogun/base/Parallel.h>
19 
25 
26 #include <vector>
27 #include <string>
28 #include <fstream>
29 #include <cmath>
30 
31 #include <assert.h>
32 
33 using namespace shogun;
34 
36  : CStringKernel<char>(0)
37 {
38  init();
40 }
41 
42 CSpectrumRBFKernel::CSpectrumRBFKernel (int32_t size, float64_t *AA_matrix_, int32_t degree_, float64_t width_)
43  : CStringKernel<char>(size), alphabet(NULL), degree(degree_), width(width_), sequences(NULL), string_features(NULL), nof_sequences(0), max_sequence_length(0)
44 {
45  init();
47 
48  target_letter_0=-1 ;
49 
51 
52  sg_memcpy(AA_matrix.matrix, AA_matrix_, 128*128*sizeof(float64_t)) ;
53 
55  SGStringList<char> string_list;
56  string_list.strings = sequences;
57  string_list.num_strings = nof_sequences;
59 
60  //string_features = new CStringFeatures<char>(sequences, nof_sequences, max_sequence_length, PROTEIN);
64 }
65 
67  CStringFeatures<char>* l, CStringFeatures<char>* r, int32_t size, float64_t* AA_matrix_, int32_t degree_, float64_t width_)
68 : CStringKernel<char>(size), alphabet(NULL), degree(degree_), width(width_), sequences(NULL), string_features(NULL), nof_sequences(0), max_sequence_length(0)
69 {
70  target_letter_0=-1 ;
71 
73  sg_memcpy(AA_matrix.matrix, AA_matrix_, 128*128*sizeof(float64_t)) ;
74 
75  init(l, r);
77 }
78 
80 {
81  cleanup();
83  SG_FREE(sequences);
84 }
85 
87 {
88 
89  int32_t aa_to_index[128];//profile
90  aa_to_index[(uint8_t) 'A'] = 0;
91  aa_to_index[(uint8_t) 'R'] = 1;
92  aa_to_index[(uint8_t) 'N'] = 2;
93  aa_to_index[(uint8_t) 'D'] = 3;
94  aa_to_index[(uint8_t) 'C'] = 4;
95  aa_to_index[(uint8_t) 'Q'] = 5;
96  aa_to_index[(uint8_t) 'E'] = 6;
97  aa_to_index[(uint8_t) 'G'] = 7;
98  aa_to_index[(uint8_t) 'H'] = 8;
99  aa_to_index[(uint8_t) 'I'] = 9;
100  aa_to_index[(uint8_t) 'L'] = 10;
101  aa_to_index[(uint8_t) 'K'] = 11;
102  aa_to_index[(uint8_t) 'M'] = 12;
103  aa_to_index[(uint8_t) 'F'] = 13;
104  aa_to_index[(uint8_t) 'P'] = 14;
105  aa_to_index[(uint8_t) 'S'] = 15;
106  aa_to_index[(uint8_t) 'T'] = 16;
107  aa_to_index[(uint8_t) 'W'] = 17;
108  aa_to_index[(uint8_t) 'Y'] = 18;
109  aa_to_index[(uint8_t) 'V'] = 19;
110  SG_DEBUG("initializing background\n")
111  double background[20]; // profile
112  background[0]=0.0799912015849807; //A
113  background[1]=0.0484482507611578;//R
114  background[2]=0.044293531582512;//N
115  background[3]=0.0578891399707563;//D
116  background[4]=0.0171846021407367;//C
117  background[5]=0.0380578923048682;//Q
118  background[6]=0.0638169929675978;//E
119  background[7]=0.0760659374742852;//G
120  background[8]=0.0223465499452473;//H
121  background[9]=0.0550905793661343;//I
122  background[10]=0.0866897071203864;//L
123  background[11]=0.060458245507428;//K
124  background[12]=0.0215379186368154;//M
125  background[13]=0.0396348024787477;//F
126  background[14]=0.0465746314476874;//P
127  background[15]=0.0630028230885602;//S
128  background[16]=0.0580394726014824;//T
129  background[17]=0.0144991866213453;//W
130  background[18]=0.03635438623143;//Y
131  background[19]=0.0700241481678408;//V
132 
133 
134  std::vector<std::string> seqs;
135  //int32_t nof_sequences = 7329;
136 
137  double C = 0.8;
138  const char *filename="/fml/ag-raetsch/home/toussaint/scp/aawd_compbio_workshop/code_nora/data/profile/profiles";
139  std::ifstream fin(filename);
140 
141  SG_DEBUG("Reading profiles from %s\n", filename)
142  std::string line;
143  while (!fin.eof())
144  {
145  std::getline(fin, line);
146 
147  if (line[0] == '>') // new sequence
148  {
149  int idx = line.find_first_of(' ');
150  sequence_labels.push_back(line.substr(1,idx-1));
151  std::getline(fin, line);
152  std::string orig_sequence = line;
153  std::string sequence="";
154 
155  int len_line = line.length();
156 
157  // skip 3 lines
158 
159  std::getline(fin, line);
160  std::getline(fin, line);
161  std::getline(fin, line);
162 
163  profiles.push_back(std::vector<double>());
164 
165  std::vector<double>& curr_profile = profiles.back();
166  for (int i=0; i < len_line; ++i)
167  {
168  std::getline(fin, line);
169  int a = line.find_first_not_of(' '); // index position
170  int b = line.find_first_of(' ', a); // index position
171  a = line.find_first_not_of(' ', b); // aa position
172  b = line.find_first_of(' ', a); // aa position
173  std::string aa=line.substr(a,b-a);
174  if (0) //(aa =="B" || aa == "X" || aa == "Z")
175  {
176  int pos = seqs.size()+1;
177  SG_DEBUG("Skipping aa in sequence %d\n", pos)
178  continue;
179  }
180  else
181  {
182  sequence += aa;
183 
184  a = line.find_first_not_of(' ', b); // beginning of block to ignore
185  b = line.find_first_of(' ', a); // aa position
186 
187  for (int j=0; j < 19; ++j)
188  {
189  a = line.find_first_not_of(' ', b);
190  b = line.find_first_of(' ', a);
191  }
192 
193  int all_zeros = 1;
194  // interesting block
195  for (int j=0; j < 20; ++j)
196  {
197  a = line.find_first_not_of(' ', b);
198  b = line.find_first_of(' ', a);
199  double p = atof(line.substr(a, b-a).c_str());
200  if (p > 0)
201  {
202  all_zeros = 0;
203  }
204  double value = -1* std::log(C*(p/100)+(1-C)*background[j]); // taken from Leslie's example, C actually corresponds to 1/(1+C)
205  curr_profile.push_back(value);
206  //SG_DEBUG("seq %d aa %d value %f p %f bg %f\n", i, j, value,p, background[j])
207  }
208 
209  if (all_zeros)
210  {
211  SG_DEBUG(">>>>>>>>>>>>>>> all zeros")
212  if (aa !="B" && aa != "X" && aa != "Z")
213  {
214  //profile[i][temp_profile_index]=-log(C+(1-C)*background[re_candidate[temp_profile_index]]);
215  int32_t aa_index = aa_to_index[(int)aa.c_str()[0]];
216  double value = -1* std::log(C+(1-C)*background[aa_index]); // taken from Leslie's example, C actually corresponds to 1/(1+C)
217  SG_DEBUG("before %f\n", profiles.back()[(i-1) * 20 + aa_index])
218  curr_profile[(i*20) + aa_index] = value;
219  SG_DEBUG(">>> aa %c \t %d \t %f\n", aa.c_str()[0], aa_index, value)
220 
221  /*
222  for (int z=0; z <20; ++z)
223  {
224  SG_DEBUG(" %d \t %f\t", z, curr_profile[z])
225  }
226  SG_DEBUG("\n")
227  */
228  }
229  }
230  }
231  }
232 
233  if (curr_profile.size() != 20 * sequence.length())
234  {
235  SG_ERROR("Something's wrong with the profile.\n")
236  break;
237  }
238 
239  seqs.push_back(sequence);
240 
241 
242  /*
243  // 6 irrelevant lines
244  for (int i=0; i < 6; ++i)
245  {
246  std::getline(fin, line);
247  }
248  //
249  */
250  }
251  }
252 
253  fin.close();
254 
255  nof_sequences = seqs.size();
256  sequences = SG_MALLOC(SGString<char>, nof_sequences);
257 
258  int max_len = 0;
259  for (int i=0; i < nof_sequences; ++i)
260  {
261  int len = seqs[i].length();
262  sequences[i].string = SG_MALLOC(char, len+1);
263  sequences[i].slen = len;
264  strcpy(sequences[i].string, seqs[i].c_str());
265 
266  if (len > max_len) max_len = len;
267  }
268 
269  max_sequence_length = max_len;
270  //string_features = new CStringFeatures<char>(sequences, nof_sequences, max_sequence_length, PROTEIN);
271 
272 }
273 
274 bool CSpectrumRBFKernel::init(CFeatures* l, CFeatures* r)
275 {
276  // >> profile
277 /*
278  read_profiles_and_sequences();
279  l = string_features;
280  r = string_features;
281  */
282  // << profile
283 
284  int32_t lhs_changed=(lhs!=l);
285  int32_t rhs_changed=(rhs!=r);
286 
288 
289  SG_DEBUG("lhs_changed: %i\n", lhs_changed)
290  SG_DEBUG("rhs_changed: %i\n", rhs_changed)
291 
294 
296  alphabet=sf_l->get_alphabet();
297  CAlphabet* ralphabet=sf_r->get_alphabet();
298 
299  if (!((alphabet->get_alphabet()==DNA) || (alphabet->get_alphabet()==RNA)))
300  properties &= ((uint64_t) (-1)) ^ (KP_LINADD | KP_BATCHEVALUATION);
301 
302  ASSERT(ralphabet->get_alphabet()==alphabet->get_alphabet())
303  SG_UNREF(ralphabet);
304 
305 
306  return init_normalizer();
307 }
308 
310 {
311 
313  alphabet=NULL;
314 
316 }
317 
318 inline bool isaa(char c)
319 {
320  if (c<65 || c>89 || c=='B' || c=='J' || c=='O' || c=='U' || c=='X' || c=='Z')
321  return false ;
322  return true ;
323 }
324 
325 float64_t CSpectrumRBFKernel::AA_helper(const char* path, const int seq_degree, const char* joint_seq, unsigned int index)
326 {
327  //const char* AA = "ARNDCQEGHILKMFPSTWYV";
328  float64_t diff=0.0 ;
329 
330  for (int i=0; i<seq_degree; i++)
331  {
332  if (!isaa(path[i])||!isaa(joint_seq[index+i]))
333  diff+=1.4 ;
334  else
335  {
336  diff += AA_matrix.matrix[ (path[i]-1)*128 + path[i] - 1] ;
337  diff -= 2*AA_matrix.matrix[ (path[i]-1)*128 + joint_seq[index+i] - 1] ;
338  diff += AA_matrix.matrix[ (joint_seq[index+i]-1)*128 + joint_seq[index+i] - 1] ;
339  if (CMath::is_nan(diff))
340  fprintf(stderr, "nan occurred: '%c' '%c'\n", path[i], joint_seq[index+i]) ;
341  }
342  }
343 
344  return exp( - diff/width) ;
345 }
346 
347 float64_t CSpectrumRBFKernel::compute(int32_t idx_a, int32_t idx_b)
348 {
349  int32_t alen, blen;
350  bool afree, bfree;
351 
352  char* avec = ((CStringFeatures<char>*) lhs)->get_feature_vector(idx_a, alen, afree);
353  char* bvec = ((CStringFeatures<char>*) rhs)->get_feature_vector(idx_b, blen, bfree);
354 
355  float64_t result=0;
356  for (int32_t i=0; i<alen; i++)
357  {
358  for (int32_t j=0; j<blen; j++)
359  {
360  if ((i+degree<=alen) && (j+degree<=blen))
361  result += AA_helper(&(avec[i]), degree, bvec, j) ;
362  }
363  }
364 
365  ((CStringFeatures<char>*) lhs)->free_feature_vector(avec, idx_a, afree);
366  ((CStringFeatures<char>*) rhs)->free_feature_vector(bvec, idx_b, bfree);
367  return result;
368 }
369 
371  float64_t* AA_matrix_)
372 {
373 
374  if (AA_matrix_)
375  {
376  SG_DEBUG("Setting AA_matrix\n")
377  sg_memcpy(AA_matrix.matrix, AA_matrix_, 128*128*sizeof(float64_t)) ;
378  return true ;
379  }
380 
381  return false;
382 }
383 
385 {
386  SG_ADD(&degree, "degree", "degree of the kernel", MS_AVAILABLE);
387  SG_ADD(&AA_matrix, "AA_matrix", "128*128 scalar product matrix", MS_NOT_AVAILABLE);
388  SG_ADD(&width, "width", "width of Gaussian", MS_AVAILABLE);
389  SG_ADD(&nof_sequences, "nof_sequences", "length of the sequence",
391  m_parameters->add_vector(&sequences, &nof_sequences, "sequences", "the sequences as a part of profile");
393  "max_sequence_length", "max length of the sequence", MS_NOT_AVAILABLE);
394 }
395 
397 {
398  SG_ADD((CSGObject**)&alphabet, "alphabet", "the alphabet used by kernel",
400 }
401 
402 void CSpectrumRBFKernel::init()
403 {
404  alphabet = NULL;
405  degree = 0;
406  width = 0.0;
407  sequences = NULL;
408  string_features = NULL;
409  nof_sequences = 0;
411 
412  initialized = false;
413 
414  max_mismatch = 0;
415  target_letter_0 = 0;
416 }
RNA - letters A,C,G,U.
Definition: Alphabet.h:32
virtual void cleanup()
Definition: Kernel.cpp:171
template class SGStringList
Definition: SGObject.h:41
DNA - letters A,C,G,T.
Definition: Alphabet.h:26
virtual bool init(CFeatures *l, CFeatures *r)
EAlphabet get_alphabet() const
Definition: Alphabet.h:130
#define SG_ERROR(...)
Definition: SGIO.h:128
The class Alphabet implements an alphabet and alphabet utility functions.
Definition: Alphabet.h:91
Parameter * m_parameters
Definition: SGObject.h:567
float64_t AA_helper(const char *path, const int degree, const char *joint_seq, unsigned int index)
CStringFeatures< char > * string_features
IUPAC_AMINO_ACID.
Definition: Alphabet.h:53
bool set_AA_matrix(float64_t *AA_matrix_)
#define SG_REF(x)
Definition: SGObject.h:52
#define ASSERT(x)
Definition: SGIO.h:200
Class SGObject is the base class of all shogun objects.
Definition: SGObject.h:125
SGString< char > * sequences
std::vector< std::vector< float64_t > > profiles
double float64_t
Definition: common.h:60
virtual bool init_normalizer()
Definition: Kernel.cpp:166
bool isaa(char c)
CFeatures * rhs
feature vectors to occur on right hand side
#define SG_UNREF(x)
Definition: SGObject.h:53
void add_vector(bool **param, index_t *length, const char *name, const char *description="")
Definition: Parameter.cpp:335
#define SG_DEBUG(...)
Definition: SGIO.h:106
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
CFeatures * lhs
feature vectors to occur on left hand side
static int is_nan(double f)
checks whether a float is nan
Definition: Math.cpp:228
The class Features is the base class of all feature objects.
Definition: Features.h:68
SGString< T > * strings
Definition: SGStringList.h:88
index_t slen
Definition: SGString.h:79
#define SG_ADD(...)
Definition: SGObject.h:94
float64_t compute(int32_t idx_a, int32_t idx_b)
Template class StringKernel, is the base class of all String Kernels.
Definition: StringKernel.h:26
SGMatrix< float64_t > AA_matrix
std::vector< std::string > sequence_labels

SHOGUN Machine Learning Toolbox - Documentation