GPUMLib  0.2.2
GPU Machine Learning Library
rbf.h
1 /*
2  Ricardo Quintas is a MSc Student at the University of Coimbra, Portugal
3  Copyright (C) 2009, 2010 Ricardo Quintas
4 
5  This file is part of GPUMLib.
6 
7  GPUMLib is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  This program is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with this program. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
23 
25 #ifndef RBF_h
26 #define RBF_h
27 
28 #include <vector>
29 #include <algorithm>
30 
31 #include "../Common/CudaDefinitions.h"
32 #include "../memory/HostMatrix.h"
33 #include "../memory/DeviceMatrix.h"
34 
35 #include "utils.h"
36 #include "RBFkernels.h"
37 
38 #define imin(X, Y) ((X) < (Y) ? (X) : (Y))
39 
40 using namespace std;
41 
42 unsigned int seed = (unsigned)time(0);
43 
44 
45 float *c_width;
46 
47 float beta = 3;
48 
49 
50 typedef std::pair<int, float> my_pair;
51 
55 void arrayShuffle(int *array, int length) {
56  for(int i = 0; i < length; i++) {
57  int r = rand() % length;
58  int t = array[r];
59  array[r] = array[i];
60  array[i] = t;
61  }
62 }
63 
65 bool sort_pred(const my_pair& left, const my_pair& right)
66 {
67  return left.second < right.second;
68 }
69 
70 
72 int randint(int lowest, int highest){
73 
74  srand(seed);
75 
76  int range=(highest-lowest)+1;
77 
78  return lowest+int(range*rand()/(RAND_MAX + 1.0));
79 }
80 
87 double euclidian_distance(const Eigen::MatrixXd &A, int idxA, const Eigen::MatrixXd &B, int idxB){
88 
89  int i;
90 
91  double sum = 0;
92 
93  double a;
94  double b;
95 
96  for(i = 0; i < B.rows(); i++){
97 
98  a = A(idxA,i);
99  b = B(i,idxB);
100 
101  sum = sum + pow( a - b , 2);
102 
103  }
104 
105  return sqrt(sum);
106 }
107 
115 float basis_function(const Eigen::MatrixXd &A, int idxA, const Eigen::MatrixXd &B, int idxB,
116  double (*f)(const Eigen::MatrixXd&, int,const Eigen::MatrixXd&, int)){
117 
118  //return exp(-c_width[idxA] * pow(f(A,idxA,B,idxB),2));
119 
120  return (float) exp(-(pow(f(A,idxA,B,idxB),2)/pow(c_width[idxA],2)));
121 
122 }
123 
128 Eigen::MatrixXd activation(const Eigen::MatrixXd &Input, const Eigen::MatrixXd &Centers){
129 
130  int i, j;
131 
132  Eigen::MatrixXd Output(Input.cols(),Centers.rows());
133 
134  for(i = 0; i < Centers.rows(); i++){
135 
136  for(j = 0; j < Input.cols() ; j++){
137 
138  Output(j,i) = basis_function(Centers,i,Input,j,euclidian_distance);
139 
140  }
141 
142  }
143 
144  return Output;
145 }
146 
152 Eigen::MatrixXd trainRBF(const Eigen::MatrixXd &Input,const Eigen::MatrixXd &Target,const Eigen::MatrixXd &Centers){
153 
154  Eigen::MatrixXd G = activation(Input,Centers);
155  Eigen::MatrixXd Aplus = UTILS::pseudoinverse(G);
156 
157  return Aplus*Target.transpose();
158 }
159 
165 Eigen::MatrixXd testRBF(const Eigen::MatrixXd &Input,const Eigen::MatrixXd &Centers, const Eigen::MatrixXd &Weights){
166  Eigen::MatrixXd G = activation(Input,Centers);
167 
168  return G*Weights;
169 }
170 
176 Eigen::MatrixXd addCenter(const Eigen::MatrixXd &X,const Eigen::MatrixXd &Centers, int index){
177 
178  Eigen::MatrixXd aux = UTILS::resize(Centers,Centers.rows()+1,Centers.cols());
179 
180  for(int i = 0; i < X.rows(); i++)
181  aux(aux.rows()-1,i) = X(i,index);
182 
183  return aux;
184 }
185 
193 void rbf(Eigen::MatrixXd Input,Eigen::MatrixXd Target,int network_size, int number_neighbours,Eigen::MatrixXd &OutWeights,Eigen::MatrixXd &OutCenters){
194 
195  int i,j;
196 
197  Eigen::MatrixXd Output(Input.rows(),1);
198 
199  c_width = (float*) malloc(sizeof(float)*network_size);
200  memset(c_width,beta,sizeof(float)*network_size);
201 
202  /*Random initialization of centers*/
203 
204  int *used;
205  used = (int*) malloc(sizeof(int)*Target.cols());
206 
207  for(i = 0; i < Target.cols(); i++) {
208  used[i]= i;
209  }
210  arrayShuffle(used, Target.cols());
211 
212  Eigen::MatrixXd Centers = Eigen::MatrixXd(1,Input.rows());
213 
214  for(int i = 0; i < Input.rows(); i++)
215  Centers(0,i) = Input(i,used[0]);
216 
217  for(i = 1; i < network_size; i++){
218  Centers = addCenter(Input,Centers,used[i]);
219  }
220 
221  free(used);
222 
223  /*Selection of centers with k-means*/
224  bool changed = true;
225  int *count = (int*) malloc(sizeof(int)*Centers.rows());
226 
227  HostMatrix<float> host_X = HostMatrix<float>(Input.rows(),Input.cols());
228 
229  int k = 0;
230  for(i = 0; i < Input.rows(); i++){
231  for(j = 0; j < Input.cols(); j++){
232  host_X.Pointer()[k] = Input(i,j);
233  k++;
234  }
235  }
236 
237  DeviceMatrix<float> device_X = DeviceMatrix<float>(host_X);
238 
239  HostMatrix<float> host_Centers = HostMatrix<float>(Centers.rows(),Centers.cols());
240 
241  k = 0;
242  for(i = 0; i < Centers.rows(); i++){
243  for(j = 0; j < Centers.cols(); j++){
244  host_Centers.Pointer()[k] = Centers(i,j);
245  k++;
246  }
247  }
248 
249  DeviceMatrix<float> device_Centers = DeviceMatrix<float>(host_Centers);
250 
251  HostMatrix<float> host_output2 = HostMatrix<float>(Input.cols(),Centers.rows());
252  DeviceMatrix<float> device_output2 = DeviceMatrix<float>(host_output2);
253 
254  int *attrib_center;
255  attrib_center = (int*) malloc(sizeof(int)*host_output2.Rows());
256 
257  int *device_attrib_center;
258  cudaMalloc((void**)&device_attrib_center, sizeof(int)*host_output2.Rows());
259 
260  while(changed){
261 
262  host_Centers = HostMatrix<float>(Centers.rows(),Centers.cols());
263 
264  k = 0;
265  for(i = 0; i < Centers.rows(); i++){
266  for(j = 0; j < Centers.cols(); j++){
267  host_Centers.Pointer()[k] = Centers(i,j);
268  k++;
269  }
270  }
271 
272  device_Centers = DeviceMatrix<float>(host_Centers);
273 
274  /*Check distances between data and Centers*/
275  KernelEuclidianDistance(device_output2.Pointer(), device_output2.Rows(), device_output2.Columns(), device_X.Pointer(), device_X.Columns(), device_Centers.Pointer(), device_Centers.Columns());
276 
277  /*Attribution of samples to centers*/
278  memset(attrib_center,0,sizeof(int)*host_output2.Rows());
279  cudaMemset(device_attrib_center,0,sizeof(int)*host_output2.Rows());
280 
281  KernelCenterAttribution(device_output2.Pointer(),device_output2.Rows(), device_output2.Columns(),device_attrib_center);
282  cudaMemcpy(attrib_center, device_attrib_center, sizeof(int)*host_output2.Rows(), cudaMemcpyDeviceToHost);
283 
284 
285  changed = false;
286  memset(count,0,sizeof(int)*host_output2.Columns());
287 
288  Eigen::MatrixXd newCenters = Eigen::MatrixXd(network_size,Input.rows());
289  newCenters.setZero();
290 
291  /*Copy data to new centers, averaging data*/
292  for(i = 0; i < host_output2.Rows(); i++){
293 
294  for(j = 0; j < newCenters.cols(); j++){
295 
296  newCenters(attrib_center[i],j) = (newCenters(attrib_center[i],j) * count[attrib_center[i]] + Input(j,i))/(count[attrib_center[i]]+1);
297  }
298 
299  count[attrib_center[i]] += 1;
300 
301  }
302 
303  /*Check if centers changed*/
304 
305 
306  for(i = 0; i < Centers.rows(); i++){
307  if(euclidian_distance(Centers,i,newCenters.transpose(),i) > 0){
308  changed = true;
309  break;
310  }
311  }
312 
313  Centers = newCenters;
314 
315  }
316 
317  free(attrib_center);
318 
319  free(count);
320 
321  /*Adjust width using mean of distance to neighbours*/
322  for(i = 0; i < Centers.rows(); i++){
323 
324  std::vector<std::pair<int,float>> tmp;
325 
326  for(j = 0; j < Centers.rows(); j++){
327  tmp.push_back(std::pair<int,float>(j,euclidian_distance(Centers,i,Centers.transpose(),j)));
328  }
329 
330  std::sort(tmp.begin(), tmp.end(), sort_pred);
331 
332  c_width[i] = 0;
333 
334  for(j=0;j<=number_neighbours;j++){
335  std::pair<int,float> aux = tmp[j];
336  c_width[i] += aux.second;
337  }
338 
339  c_width[i] = c_width[i]/number_neighbours;
340  }
341 
342  /*Training*/
343  Eigen::MatrixXd Weights = trainRBF(Input,Target,Centers);
344 
345 
346  /*Return Weights and Centers*/
347  OutWeights = Weights;
348  OutCenters = Centers;
349 
350 }
351 
352 #endif
353 
354 //! @}
Eigen::MatrixXd activation(const Eigen::MatrixXd &Input, const Eigen::MatrixXd &Centers)
Definition: rbf.h:128
void arrayShuffle(int *array, int length)
Definition: rbf.h:55
float basis_function(const Eigen::MatrixXd &A, int idxA, const Eigen::MatrixXd &B, int idxB, double(*f)(const Eigen::MatrixXd &, int, const Eigen::MatrixXd &, int))
Definition: rbf.h:115
int randint(int lowest, int highest)
Generates a random Integer between lowest and highest values.
Definition: rbf.h:72
Eigen::MatrixXd trainRBF(const Eigen::MatrixXd &Input, const Eigen::MatrixXd &Target, const Eigen::MatrixXd &Centers)
Definition: rbf.h:152
void rbf(Eigen::MatrixXd Input, Eigen::MatrixXd Target, int network_size, int number_neighbours, Eigen::MatrixXd &OutWeights, Eigen::MatrixXd &OutCenters)
Definition: rbf.h:193
double euclidian_distance(const Eigen::MatrixXd &A, int idxA, const Eigen::MatrixXd &B, int idxB)
Definition: rbf.h:87
Type * Pointer() const
Definition: BaseMatrix.h:88
int Columns() const
Definition: BaseMatrix.h:80
void KernelCenterAttribution(cudafloat *Output, int output_height, int output_width, int *attrib_center)
Eigen::MatrixXd addCenter(const Eigen::MatrixXd &X, const Eigen::MatrixXd &Centers, int index)
Definition: rbf.h:176
bool sort_pred(const my_pair &left, const my_pair &right)
Sorting predicate, for use with the sort function.
Definition: rbf.h:65
int Rows() const
Definition: BaseMatrix.h:74
void KernelEuclidianDistance(cudafloat *d_C, cudafloat *d_A, cudafloat *d_B, int uiWA, int uiWB, int uiWC, int uiHC)
Eigen::MatrixXd testRBF(const Eigen::MatrixXd &Input, const Eigen::MatrixXd &Centers, const Eigen::MatrixXd &Weights)
Definition: rbf.h:165