SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SNPFeatures.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) 2009 Soeren Sonnenburg
8  * Copyright (C) 2009 Fraunhofer Institute FIRST and Max-Planck-Society
9  */
10 
12 #include <shogun/io/SGIO.h>
14 #include <shogun/lib/memory.h>
15 
16 using namespace shogun;
17 
19 {
20  SG_UNSTABLE("CSNPFeatures::CSNPFeatures()", "\n")
21 
22  strings = NULL;
23 
24  string_length = 0;
25  num_strings = 0;
26  w_dim = 0;
27 
28  normalization_const = 0.0;
29 
30  m_str_min = NULL;
31  m_str_maj = NULL;
32 }
33 
35  m_str_min(NULL), m_str_maj(NULL)
36 {
37  ASSERT(str)
38  ASSERT(str->have_same_length())
39  SG_REF(str);
40 
41  strings=str;
43  ASSERT((string_length & 1) == 0) // length divisible by 2
44  w_dim=3*string_length/2;
46  CAlphabet* alpha=str->get_alphabet();
47  ASSERT(alpha->get_alphabet()==SNP)
48  SG_UNREF(alpha);
49 
52 
53 }
54 
56  : CDotFeatures(orig), strings(orig.strings),
57  normalization_const(orig.normalization_const),
58  m_str_min(NULL), m_str_maj(NULL)
59 {
60  SG_REF(strings);
61 
62  if (strings)
63  {
65  ASSERT((string_length & 1) == 0) // length divisible by 2
68  }
69  else
70  {
71  string_length = 0;
72  w_dim = 0;
73  num_strings = 0;
74  }
75 
77 }
78 
80 {
82 }
83 
85 {
86  return w_dim;
87 }
88 
90 {
91  return w_dim/3;
92 }
93 
95 {
96  return F_UNKNOWN;
97 }
98 
100 {
101  return C_WD;
102 }
103 
105 {
106  return num_strings;
107 }
108 
110 {
111  return normalization_const;
112 }
113 
115 {
116  m_str_min=(uint8_t*) get_strdup(str);
117 }
118 
120 {
121  m_str_maj=(uint8_t*) get_strdup(str);
122 }
123 
125 {
126  return (char*) m_str_min;
127 }
128 
130 {
131  return (char*) m_str_maj;
132 }
133 
134 float64_t CSNPFeatures::dot(int32_t idx_a, CDotFeatures* df, int32_t idx_b)
135 {
136  ASSERT(df)
139  CSNPFeatures* sf=(CSNPFeatures*) df;
140 
141  int32_t alen, blen;
142  bool free_avec, free_bvec;
143 
144  uint8_t* avec = strings->get_feature_vector(idx_a, alen, free_avec);
145  uint8_t* bvec = sf->strings->get_feature_vector(idx_b, blen, free_bvec);
146 
147  ASSERT(alen==blen)
148  if (alen!=string_length)
149  SG_ERROR("alen (%d) !=string_length (%d)\n", alen, string_length)
152 
153  float64_t total=0;
154 
155  for (int32_t i = 0; i<alen-1; i+=2)
156  {
157  int32_t sumaa=0;
158  int32_t sumbb=0;
159  int32_t sumab=0;
160 
161  uint8_t a1=avec[i];
162  uint8_t a2=avec[i+1];
163  uint8_t b1=bvec[i];
164  uint8_t b2=bvec[i+1];
165 
166  if ((a1!=a2 || a1=='0' || a1=='0') && (b1!=b2 || b1=='0' || b2=='0'))
167  sumab++;
168  else if (a1==a2 && b1==b2)
169  {
170  if (a1!=b1)
171  continue;
172 
173  if (a1==m_str_min[i])
174  sumaa++;
175  else if (a1==m_str_maj[i])
176  sumbb++;
177  else
178  {
179  SG_ERROR("The impossible happened i=%d a1=%c "
180  "a2=%c b1=%c b2=%c min=%c maj=%c\n", i, a1,a2, b1,b2, m_str_min[i], m_str_maj[i]);
181  }
182 
183  }
184  total+=sumaa+sumbb+sumab;
185  }
186 
187  strings->free_feature_vector(avec, idx_a, free_avec);
188  sf->strings->free_feature_vector(bvec, idx_b, free_bvec);
189  return total;
190 }
191 
192 float64_t CSNPFeatures::dense_dot(int32_t vec_idx1, const float64_t* vec2, int32_t vec2_len)
193 {
194  if (vec2_len != w_dim)
195  SG_ERROR("Dimensions don't match, vec2_dim=%d, w_dim=%d\n", vec2_len, w_dim)
196 
197  float64_t sum=0;
198  int32_t len;
199  bool free_vec1;
200  uint8_t* vec = strings->get_feature_vector(vec_idx1, len, free_vec1);
201  int32_t offs=0;
202 
203  for (int32_t i=0; i<len; i+=2)
204  {
205  int32_t dim=0;
206 
207  char a1=vec[i];
208  char a2=vec[i+1];
209 
210  if (a1==a2 && a1!='0' && a2!='0')
211  {
212  if (a1==m_str_min[i])
213  dim=1;
214  else if (a1==m_str_maj[i])
215  dim=2;
216  else
217  {
218  SG_ERROR("The impossible happened i=%d a1=%c a2=%c min=%c maj=%c\n",
219  i, a1,a2, m_str_min[i], m_str_maj[i]);
220  }
221  }
222 
223  sum+=vec2[offs+dim];
224  offs+=3;
225  }
226  strings->free_feature_vector(vec, vec_idx1, free_vec1);
227 
228  return sum/normalization_const;
229 }
230 
231 void CSNPFeatures::add_to_dense_vec(float64_t alpha, int32_t vec_idx1, float64_t* vec2, int32_t vec2_len, bool abs_val)
232 {
233  if (vec2_len != w_dim)
234  SG_ERROR("Dimensions don't match, vec2_dim=%d, w_dim=%d\n", vec2_len, w_dim)
235 
236  int32_t len;
237  bool free_vec1;
238  uint8_t* vec = strings->get_feature_vector(vec_idx1, len, free_vec1);
239  int32_t offs=0;
240 
241  if (abs_val)
242  alpha=CMath::abs(alpha);
243 
244  for (int32_t i=0; i<len; i+=2)
245  {
246  int32_t dim=0;
247 
248  char a1=vec[i];
249  char a2=vec[i+1];
250 
251  if (a1==a2 && a1!='0' && a2!='0')
252  {
253  if (a1==m_str_min[i])
254  dim=1;
255  else if (a1==m_str_maj[i])
256  dim=2;
257  else
258  {
259  SG_ERROR("The impossible happened i=%d a1=%c a2=%c min=%c maj=%c\n",
260  i, a1,a2, m_str_min[i], m_str_maj[i]);
261  }
262  }
263 
264  vec2[offs+dim]+=alpha;
265  offs+=3;
266  }
267  strings->free_feature_vector(vec, vec_idx1, free_vec1);
268 }
269 
270 void CSNPFeatures::find_minor_major_strings(uint8_t* minor, uint8_t* major)
271 {
272  for (int32_t i=0; i<num_strings; i++)
273  {
274  int32_t len;
275  bool free_vec;
276  uint8_t* vec = ((CStringFeatures<uint8_t>*) strings)->get_feature_vector(i, len, free_vec);
277  ASSERT(string_length==len)
278 
279  for (int32_t j=0; j<len; j++)
280  {
281  // skip sequencing errors
282  if (vec[j]=='0')
283  continue;
284 
285  if (minor[j]==0)
286  minor[j]=vec[j];
287  else if (major[j]==0 && vec[j]!=minor[j])
288  major[j]=vec[j];
289  }
290 
291  ((CStringFeatures<uint8_t>*) strings)->free_feature_vector(vec, i, free_vec);
292  }
293 }
294 
296 {
297  SG_FREE(m_str_min);
298  SG_FREE(m_str_maj);
299  size_t tlen=(string_length+1)*sizeof(uint8_t);
300 
301  m_str_min=SG_CALLOC(uint8_t, tlen);
302  m_str_maj=SG_CALLOC(uint8_t, tlen);
303 
304  find_minor_major_strings(m_str_min, m_str_maj);
305 
306  if (snp)
307  snp->find_minor_major_strings(m_str_min, m_str_maj);
308 
309  for (int32_t j=0; j<string_length; j++)
310  {
311  // if only one symbol occurs use 0
312  if (m_str_min[j]==0)
313  m_str_min[j]='0';
314  if (m_str_maj[j]==0)
315  m_str_maj[j]='0';
316 
317  if (m_str_min[j]>m_str_maj[j])
319  }
320 }
321 
323 {
324  if (n==0)
325  {
328  }
329  else
331 
332  SG_DEBUG("normalization_const:%f\n", normalization_const)
333 }
334 
335 void* CSNPFeatures::get_feature_iterator(int32_t vector_index)
336 {
337  return NULL;
338 }
339 
340 bool CSNPFeatures::get_next_feature(int32_t& index, float64_t& value, void* iterator)
341 {
342  return false;
343 }
344 
346 {
347 }
348 
350 {
351  return new CSNPFeatures(*this);
352 }
353 
355 {
356  int32_t nsym=3;
357  float64_t* h= SG_CALLOC(float64_t, size_t(nsym)*string_length/2);
358 
359  float64_t* h_normalizer=SG_MALLOC(float64_t, string_length/2);
360  memset(h_normalizer, 0, string_length/2*sizeof(float64_t));
361  int32_t num_str=get_num_vectors();
362  for (int32_t i=0; i<num_str; i++)
363  {
364  int32_t len;
365  bool free_vec;
366  uint8_t* vec = strings->get_feature_vector(i, len, free_vec);
367 
368  for (int32_t j=0; j<len; j+=2)
369  {
370  int32_t dim=0;
371 
372  char a1=vec[j];
373  char a2=vec[j+1];
374 
375  if (a1==a2 && a1!='0' && a2!='0')
376  {
377  if (a1==m_str_min[j])
378  dim=1;
379  else if (a1==m_str_maj[j])
380  dim=2;
381  else
382  {
383  SG_ERROR("The impossible happened j=%d a1=%c a2=%c min=%c maj=%c\n",
384  j, a1,a2, m_str_min[j], m_str_maj[j]);
385  }
386  }
387 
388  h[int64_t(j/2)*nsym+dim]++;
389  h_normalizer[j/2]++;
390  }
391 
392  strings->free_feature_vector(vec, i, free_vec);
393  }
394 
395  if (normalize)
396  {
397  for (int32_t i=0; i<string_length/2; i++)
398  {
399  for (int32_t j=0; j<nsym; j++)
400  {
401  if (h_normalizer && h_normalizer[i])
402  h[int64_t(i)*nsym+j]/=h_normalizer[i];
403  }
404  }
405  }
406  SG_FREE(h_normalizer);
407 
408  return SGMatrix<float64_t>(h, nsym, string_length/2);
409 }
410 
412 {
413 
415  int32_t len=pos->strings->get_max_vector_length();
416 
417  float64_t* table=SG_MALLOC(float64_t, 3*2*len/2);
418 
419  SGMatrix<float64_t> p_hist=pos->get_histogram(false);
420  SGMatrix<float64_t> n_hist=neg->get_histogram(false);
421 
422 
423  for (int32_t i=0; i<3*len/2; i++)
424  {
425  table[2*i]=p_hist.matrix[i];
426  table[2*i+1]=n_hist.matrix[i];
427  }
428  return SGMatrix<float64_t>(table, 2,3*len/2);
429 }

SHOGUN Machine Learning Toolbox - Documentation