LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
keras_model.cc
Go to the documentation of this file.
1 #include "keras_model.h"
2 
3 #include "tbb/parallel_for.h"
4 #include "tbb/tbb.h"
5 
6 #include <iostream>
7 #include <fstream>
8 #include <chrono>
9 #include <algorithm>
10 #include <math.h>
11 using namespace std;
12 
13 
14 std::vector<float> keras::read_1d_array(std::ifstream &fin, int cols) {
15  vector<float> arr;
16  float tmp_float;
17  char tmp_char;
18  fin >> tmp_char; // for '['
19  for(int n = 0; n < cols; ++n) {
20  fin >> tmp_float;
21  arr.push_back(tmp_float);
22  }
23  fin >> tmp_char; // for ']'
24  return arr;
25 }
26 
27 void keras::DataChunk2D::read_from_file(const std::string &fname) {
28  ifstream fin(fname.c_str());
29  fin >> m_depth >> m_rows >> m_cols;
30 
31  for(int d = 0; d < m_depth; ++d) {
32  vector<vector<float> > tmp_single_depth;
33  for(int r = 0; r < m_rows; ++r) {
34  vector<float> tmp_row = keras::read_1d_array(fin, m_cols);
35  tmp_single_depth.push_back(tmp_row);
36  }
37  data.push_back(tmp_single_depth);
38  }
39  fin.close();
40 }
41 
42 
43 void keras::LayerConv2D::load_weights(std::ifstream &fin) {
44  char tmp_char = ' ';
45  string tmp_str = "";
46  float tmp_float;
47  bool skip = false;
48  fin >> m_kernels_cnt >> m_depth >> m_rows >> m_cols >> m_border_mode;
49  if (m_border_mode == "[") { m_border_mode = "valid"; skip = true; }
50 
51  cout << "LayerConv2D " << m_kernels_cnt
52  << "x" << m_depth << "x" << m_rows << "x" << m_cols << " border_mode " << m_border_mode << endl;
53  // reading kernel weights
54  for(int k = 0; k < m_kernels_cnt; ++k) {
55  vector<vector<vector<float> > > tmp_depths;
56  for(int d = 0; d < m_depth; ++d) {
57  vector<vector<float> > tmp_single_depth;
58  for(int r = 0; r < m_rows; ++r) {
59  if (!skip) { fin >> tmp_char; } // for '['
60  else { skip = false; }
61  vector<float> tmp_row;
62  for(int c = 0; c < m_cols; ++c) {
63  fin >> tmp_float;
64  //cout << tmp_float << " ";
65  tmp_row.push_back(tmp_float);
66  }
67  fin >> tmp_char; // for ']'
68  tmp_single_depth.push_back(tmp_row);
69  }
70  tmp_depths.push_back(tmp_single_depth);
71  }
72  m_kernels.push_back(tmp_depths);
73  }
74  // reading kernel biases
75  fin >> tmp_char; // for '['
76  for(int k = 0; k < m_kernels_cnt; ++k) {
77  fin >> tmp_float;
78  m_bias.push_back(tmp_float);
79  }
80  fin >> tmp_char; // for ']'
81 
82 }
83 
85  fin >> m_activation_type;
86  cout << "Activation type " << m_activation_type << endl;
87 }
88 
90  fin >> m_pool_x >> m_pool_y;
91  cout << "MaxPooling " << m_pool_x << "x" << m_pool_y << endl;
92 }
93 
94 void keras::LayerDense::load_weights(std::ifstream &fin) {
95  fin >> m_input_cnt >> m_neurons;
96  float tmp_float;
97  char tmp_char = ' ';
98  for(int i = 0; i < m_input_cnt; ++i) {
99  vector<float> tmp_n;
100  fin >> tmp_char; // for '['
101  for(int n = 0; n < m_neurons; ++n) {
102  fin >> tmp_float;
103  tmp_n.push_back(tmp_float);
104  }
105  fin >> tmp_char; // for ']'
106  m_weights.push_back(tmp_n);
107  }
108  cout << "weights " << m_weights.size() << endl;
109  fin >> tmp_char; // for '['
110  for(int n = 0; n < m_neurons; ++n) {
111  fin >> tmp_float;
112  m_bias.push_back(tmp_float);
113  }
114  fin >> tmp_char; // for ']'
115  cout << "bias " << m_bias.size() << endl;
116 
117 }
118 
119 keras::KerasModel::KerasModel(const string &input_fname) {
120  load_weights(input_fname);
121 }
122 
123 
125  vector<vector<vector<float> > > im = dc->get_3d();
126 
127  size_t csize = im[0].size();
128  size_t rsize = im[0][0].size();
129  size_t size = im.size() * csize * rsize;
130  keras::DataChunkFlat *out = new DataChunkFlat(size);
131  float * y_ret = out->get_1d_rw().data();
132  for(size_t i = 0, dst = 0; i < im.size(); ++i) {
133  for(size_t j = 0; j < csize; ++j) {
134  float * row = im[i][j].data();
135  for(size_t k = 0; k < rsize; ++k) {
136  y_ret[dst++] = row[k];
137  }
138  }
139  }
140 
141  return out;
142 }
143 
144 
146  vector<vector<vector<float> > > im = dc->get_3d();
147  vector<vector<vector<float> > > y_ret;
148  for(unsigned int i = 0; i < im.size(); ++i) {
149  vector<vector<float> > tmp_y;
150  for(unsigned int j = 0; j < (unsigned int)(im[0].size()/m_pool_x); ++j) {
151  tmp_y.push_back(vector<float>((int)(im[0][0].size()/m_pool_y), 0.0));
152  }
153  y_ret.push_back(tmp_y);
154  }
155  for(unsigned int d = 0; d < y_ret.size(); ++d) {
156  for(unsigned int x = 0; x < y_ret[0].size(); ++x) {
157  unsigned int start_x = x*m_pool_x;
158  unsigned int end_x = start_x + m_pool_x;
159  for(unsigned int y = 0; y < y_ret[0][0].size(); ++y) {
160  unsigned int start_y = y*m_pool_y;
161  unsigned int end_y = start_y + m_pool_y;
162 
163  vector<float> values;
164  for(unsigned int i = start_x; i < end_x; ++i) {
165  for(unsigned int j = start_y; j < end_y; ++j) {
166  values.push_back(im[d][i][j]);
167  }
168  }
169  y_ret[d][x][y] = *max_element(values.begin(), values.end());
170  }
171  }
172  }
174  out->set_data(y_ret);
175  return out;
176 }
177 
178 void keras::missing_activation_impl(const string &act) {
179  cout << "Activation " << act << " not defined!" << endl;
180  cout << "Please add its implementation before use." << endl;
181  exit(1);
182 }
183 
185 
186  if (dc->get_data_dim() == 3) {
187  vector<vector<vector<float> > > y = dc->get_3d();
188  if(m_activation_type == "relu") {
189  for(unsigned int i = 0; i < y.size(); ++i) {
190  for(unsigned int j = 0; j < y[0].size(); ++j) {
191  for(unsigned int k = 0; k < y[0][0].size(); ++k) {
192  if(y[i][j][k] < 0) y[i][j][k] = 0;
193  }
194  }
195  }
196  } else if(m_activation_type == "tanh") {
197  for(unsigned int i = 0; i < y.size(); ++i) {
198  for(unsigned int j = 0; j < y[0].size(); ++j) {
199  for(unsigned int k = 0; k < y[0][0].size(); ++k) {
200  y[i][j][k] = tanh(y[i][j][k]);
201  }
202  }
203  }
204  } else {
205  keras::missing_activation_impl(m_activation_type);
206  }
207 
209  out->set_data(y);
210  return out;
211 
212  } else if (dc->get_data_dim() == 1) { // flat data, use 1D
213  vector<float> y = dc->get_1d();
214  if(m_activation_type == "relu") {
215  for(unsigned int k = 0; k < y.size(); ++k) {
216  if(y[k] < 0) y[k] = 0;
217  }
218  } else if(m_activation_type == "softmax") {
219  float sum = 0.0;
220  for(unsigned int k = 0; k < y.size(); ++k) {
221  y[k] = exp(y[k]);
222  sum += y[k];
223  }
224  for(unsigned int k = 0; k < y.size(); ++k) {
225  y[k] /= sum;
226  }
227  } else if(m_activation_type == "tanh") {
228  for(unsigned int k = 0; k < y.size(); ++k) {
229  y[k] = tanh(y[k]);
230  }
231  } else if(m_activation_type == "sigmoid") {
232  for(unsigned int k = 0; k < y.size(); ++k) {
233  y[k] = 1.0F / (1.0F + exp(-y[k]));
234  }
235  } else {
236  keras::missing_activation_impl(m_activation_type);
237  }
238 
239  keras::DataChunk *out = new DataChunkFlat();
240  out->set_data(y);
241  return out;
242 
243  } else { throw "data dim not supported"; }
244 
245  return dc;
246 }
247 
248 // with border mode = valid
250  std::vector< std::vector<float> > & y, // accumulate here
251  std::vector< std::vector<float> > const & im,
252  std::vector< std::vector<float> > const & k)
253 {
254  size_t k1_size = k.size(), k2_size = k[0].size();
255  unsigned int st_x = (k1_size - 1) >> 1;
256  unsigned int st_y = (k2_size - 1) >> 1;
257 
258  for(unsigned int i = st_x; i < im.size()-st_x; ++i) {
259  for(unsigned int j = st_y; j < im[0].size()-st_y; ++j) {
260 
261  float sum = 0;
262  for(unsigned int k1 = 0; k1 < k1_size; ++k1) {
263  //const float * k_data = k[k1_size-k1-1].data();
264  //const float * im_data = im[i-st_x+k1].data();
265  for(unsigned int k2 = 0; k2 < k2_size; ++k2) {
266  //sum += k_data[k2_size-k2-1] * im_data[j-st_y+k2];
267  sum += k[k1_size-k1-1][k2_size-k2-1] * im[i-st_x+k1][j-st_y+k2];
268  }
269  }
270  y[i-st_x][j-st_y] += sum;
271  }
272  }
273 }
274 
275 // with border mode = same
277  std::vector< std::vector<float> > & y, // accumulate here
278  std::vector< std::vector<float> > const & im,
279  std::vector< std::vector<float> > const & k)
280 {
281  size_t k1_size = k.size(), k2_size = k[0].size();
282  unsigned int st_x = (k1_size - 1) >> 1;
283  unsigned int st_y = (k2_size - 1) >> 1;
284 
285  size_t max_imc = im.size() - 1;
286  size_t max_imr = im[0].size() - 1;
287 
288  for(unsigned int i = 0; i < im.size(); ++i) {
289  for(unsigned int j = 0; j < im[0].size(); ++j) {
290 
291  float sum = 0;
292  for(unsigned int k1 = 0; k1 < k1_size; ++k1) {
293  //const float * k_data = k[k1_size-k1-1].data();
294  //const float * im_data = im[i-st_x+k1].data();
295  for(unsigned int k2 = 0; k2 < k2_size; ++k2) {
296  if(i+k1 < st_x) continue;
297  if(i+k1 > max_imc+st_x) continue;
298  if(j+k2 < st_y) continue;
299  if(j+k2 > max_imr+st_y) continue;
300 
301  //sum += k_data[k2_size-k2-1] * im_data[j-st_y+k2];
302  sum += k[k1_size-k1-1][k2_size-k2-1] * im[i-st_x+k1][j-st_y+k2];
303  }
304  }
305  y[i][j] += sum;
306  }
307  }
308 }
309 
311  unsigned int st_x = (m_kernels[0][0].size()-1) >> 1;
312  unsigned int st_y = (m_kernels[0][0][0].size()-1) >> 1;
313  auto const & im = dc->get_3d();
314 
315  size_t size_x = (m_border_mode == "valid")? im[0].size() - 2 * st_x : im[0].size();
316  size_t size_y = (m_border_mode == "valid")? im[0][0].size() - 2 * st_y: im[0][0].size();
317 
318  // depth rows cols
319  keras::DataChunk2D *out = new keras::DataChunk2D(m_kernels.size(), size_x, size_y, 0);
320  auto & y_ret = out->get_3d_rw();
321 
322  //auto t1 = std::chrono::high_resolution_clock::now();
323 
324  // Parallelize the kernal calculation
325  tbb::parallel_for( size_t(0), size_t(m_kernels.size()), [&]( size_t j ) {
326 
327  for(unsigned int m = 0; m < im.size(); ++m) { // loop over image depth
328  if (m_border_mode == "valid")
329  keras::conv_single_depth_valid(y_ret[j], im[m], m_kernels[j][m]);
330  else
331  keras::conv_single_depth_same(y_ret[j], im[m], m_kernels[j][m]);
332  }
333 
334  for(unsigned int x = 0; x < y_ret[0].size(); ++x) {
335 
336  size_t size = y_ret[0][0].size();
337  float bias = m_bias[j];
338  size_t k = 0;
339 
340  for(unsigned int y = 0; y < size/8; ++y) {
341  y_ret[j][x][k] += bias;
342  y_ret[j][x][k+1] += bias;
343  y_ret[j][x][k+2] += bias;
344  y_ret[j][x][k+3] += bias;
345  y_ret[j][x][k+4] += bias;
346  y_ret[j][x][k+5] += bias;
347  y_ret[j][x][k+6] += bias;
348  y_ret[j][x][k+7] += bias;
349  k += 8;
350  }
351  while (k < size) { y_ret[j][x][k] += bias; ++k; }
352 
353  }
354  });
355 
356  //auto t2 = std::chrono::high_resolution_clock::now();
357  /*
358  cout << "Parallal : "
359  << std::chrono::duration_cast<std::chrono::microseconds>(t2-t1).count()
360  << " microseconds." << endl;
361  */
362 
363  return out;
364 }
365 
367  //cout << "weights: input size " << m_weights.size() << endl;
368  //cout << "weights: neurons size " << m_weights[0].size() << endl;
369  //cout << "bias " << m_bias.size() << endl;
370  size_t size = m_weights[0].size();
371  size_t size8 = size >> 3;
372 
373  keras::DataChunkFlat *out = new DataChunkFlat(size, 0);
374  float * y_ret = out->get_1d_rw().data();
375 
376  auto const & im = dc->get_1d();
377 
378  for (size_t j = 0; j < m_weights.size(); ++j) { // iter over input
379  const float * w = m_weights[j].data();
380  float p = im[j];
381  size_t k = 0;
382  for (size_t i = 0; i < size8; ++i) { // iter over neurons
383  y_ret[k] += w[k] * p; // vectorize if you can
384  y_ret[k+1] += w[k+1] * p;
385  y_ret[k+2] += w[k+2] * p;
386  y_ret[k+3] += w[k+3] * p;
387  y_ret[k+4] += w[k+4] * p;
388  y_ret[k+5] += w[k+5] * p;
389  y_ret[k+6] += w[k+6] * p;
390  y_ret[k+7] += w[k+7] * p;
391  k += 8;
392  }
393  while (k < size) { y_ret[k] += w[k] * p; ++k; }
394  }
395  for (size_t i = 0; i < size; ++i) { // add biases
396  y_ret[i] += m_bias[i];
397  }
398 
399  return out;
400 }
401 
402 
404  //cout << endl << "KerasModel compute output" << endl;
405  //cout << "Input data size:" << endl;
406  //dc->show_name();
407 
408  keras::DataChunk *inp = dc;
409  keras::DataChunk *out = 0;
410  for(int l = 0; l < (int)m_layers.size(); ++l) {
411  //cout << "Processing layer " << m_layers[l]->get_name() << endl;
412  out = m_layers[l]->compute_output(inp);
413 
414  //cout << "Input" << endl;
415  //inp->show_name();
416  //cout << "Output" << endl;
417  //out->show_name();
418 
419  if (inp != dc) delete inp;
420  inp = 0L;
421  inp = out;
422  }
423 
424  //cout << "Output: ";
425  //out->show_values();
426 
427  std::vector<float> flat_out = out->get_1d();
428  delete out;
429 
430  return flat_out;
431 }
432 
433 void keras::KerasModel::load_weights(const string &input_fname) {
434  cout << "Reading model from " << input_fname << endl;
435  ifstream fin(input_fname.c_str());
436  string layer_type = "";
437  string tmp_str = "";
438  int tmp_int = 0;
439 
440  fin >> tmp_str >> m_layers_cnt;
441  cout << "Layers " << m_layers_cnt << endl;
442 
443  for(int layer = 0; layer < m_layers_cnt; ++layer) { // iterate over layers
444  fin >> tmp_str >> tmp_int >> layer_type;
445  cout << "Layer " << tmp_int << " " << layer_type << endl;
446 
447  Layer *l = 0L;
448  if(layer_type == "Convolution2D") {
449  l = new LayerConv2D();
450  } else if(layer_type == "Activation") {
451  l = new LayerActivation();
452  } else if(layer_type == "MaxPooling2D") {
453  l = new LayerMaxPooling();
454  } else if(layer_type == "Flatten") {
455  l = new LayerFlatten();
456  } else if(layer_type == "Dense") {
457  l = new LayerDense();
458  } else if(layer_type == "Dropout") {
459  continue; // we dont need dropout layer in prediciton mode
460  }
461  if(l == 0L) {
462  cout << "Layer is empty, maybe it is not defined? Cannot define network." << endl;
463  return;
464  }
465  l->load_weights(fin);
466  m_layers.push_back(l);
467  }
468 
469  fin.close();
470 }
471 
473  for(int i = 0; i < (int)m_layers.size(); ++i) {
474  delete m_layers[i];
475  }
476 }
477 
479 {
480  int i = m_layers.size() - 1;
481  while ((i > 0) && (m_layers[i]->get_output_units() == 0)) --i;
482  return m_layers[i]->get_output_units();
483 }
Float_t x
Definition: compare.C:6
void load_weights(std::ifstream &fin)
Definition: keras_model.cc:43
void conv_single_depth_valid(std::vector< std::vector< float > > &y, std::vector< std::vector< float > > const &im, std::vector< std::vector< float > > const &k)
Definition: keras_model.cc:249
KerasModel(const std::string &input_fname)
Definition: keras_model.cc:119
keras::DataChunk * compute_output(keras::DataChunk *)
Definition: keras_model.cc:310
Float_t y
Definition: compare.C:6
int get_output_length() const
Definition: keras_model.cc:478
virtual void set_data(std::vector< std::vector< std::vector< float > > > const &)
Definition: keras_model.h:47
STL namespace.
void load_weights(const std::string &input_fname)
Definition: keras_model.cc:433
virtual void load_weights(std::ifstream &fin)=0
virtual std::vector< std::vector< std::vector< float > > > const & get_3d() const
Definition: keras_model.h:46
virtual size_t get_data_dim(void) const
Definition: keras_model.h:44
keras::DataChunk * compute_output(keras::DataChunk *)
Definition: keras_model.cc:366
keras::DataChunk * compute_output(keras::DataChunk *)
Definition: keras_model.cc:124
std::vector< float > compute_output(keras::DataChunk *dc)
Definition: keras_model.cc:403
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
void read_from_file(const std::string &fname)
Definition: keras_model.cc:27
TString fin
Definition: Style.C:24
std::vector< std::vector< std::vector< float > > > & get_3d_rw()
Definition: keras_model.h:62
Float_t d
Definition: plot.C:237
keras::DataChunk * compute_output(keras::DataChunk *)
Definition: keras_model.cc:184
void missing_activation_impl(const std::string &act)
Definition: keras_model.cc:178
void load_weights(std::ifstream &fin)
Definition: keras_model.cc:94
keras::DataChunk * compute_output(keras::DataChunk *)
Definition: keras_model.cc:145
std::vector< float > & get_1d_rw()
Definition: keras_model.h:102
Char_t n[5]
void load_weights(std::ifstream &fin)
Definition: keras_model.cc:84
void load_weights(std::ifstream &fin)
Definition: keras_model.cc:89
Float_t w
Definition: plot.C:23
virtual std::vector< float > const & get_1d() const
Definition: keras_model.h:45
void conv_single_depth_same(std::vector< std::vector< float > > &y, std::vector< std::vector< float > > const &im, std::vector< std::vector< float > > const &k)
Definition: keras_model.cc:276
std::vector< float > read_1d_array(std::ifstream &fin, int cols)
Definition: keras_model.cc:14