Merge branch 'tweaks' into devel
[mldemos:baraks-mldemos.git] / _AlgorithmsPlugins / KernelMethods / regressorSVR.cpp
1 /*********************************************************************\r
2 MLDemos: A User-Friendly visualization toolkit for machine learning\r
3 Copyright (C) 2010  Basilio Noris\r
4 Contact: mldemos@b4silio.com\r
5 \r
6 This library is free software; you can redistribute it and/or\r
7 modify it under the terms of the GNU Lesser General Public\r
8 License as published by the Free Software Foundation; either\r
9 version 2.1 of the License, or (at your option) any later version.\r
10 \r
11 This library is distributed in the hope that it will be useful,\r
12 but WITHOUT ANY WARRANTY; without even the implied warranty of\r
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU\r
14 Library General Public License for more details.\r
15 \r
16 You should have received a copy of the GNU Lesser General Public\r
17 License along with this library; if not, write to the Free\r
18 Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.\r
19 *********************************************************************/\r
20 #include <public.h>\r
21 #include "regressorSVR.h"\r
22 #include <nlopt/nlopt.hpp>\r
23 #include <QDebug>\r
24 \r
25 using namespace std;\r
26 \r
27 const char *RegressorSVR::GetInfoString()\r
28 {\r
29     if(!svm) return NULL;\r
30     char *text = new char[255];\r
31     sprintf(text, "%s\n", param.svm_type == NU_SVR ? "nu-SVR" : "eps-SVR");\r
32     sprintf(text, "%sKernel: ", text);\r
33     switch(param.kernel_type)\r
34     {\r
35     case LINEAR:\r
36         sprintf(text, "%s linear\n", text);\r
37         break;\r
38     case POLY:\r
39         sprintf(text, "%s polynomial (deg: %d bias: %f width: %f)\n", text, param.degree, param.coef0, param.gamma);\r
40         break;\r
41     case RBF:\r
42         sprintf(text, "%s rbf (gamma: %f)\n", text, param.gamma);\r
43         break;\r
44     case SIGMOID:\r
45         sprintf(text, "%s sigmoid (%f %f)\n", text, param.gamma, param.coef0);\r
46         break;\r
47     }\r
48     sprintf(text, "%seps: %f \t nu: %f\n", text, param.eps, param.nu);\r
49     sprintf(text, "%sSupport Vectors: %d\n", text, svm->l);\r
50     return text;\r
51 }\r
52 \r
53 RegressorSVR::RegressorSVR()\r
54     : svm(0), node(0)\r
55 {\r
56     type = REGR_SVR;\r
57     // default values\r
58     param.svm_type = EPSILON_SVR;\r
59     //param.svm_type = NU_SVR;\r
60     param.kernel_type = RBF;\r
61     param.gamma = 0.1;\r
62     param.C = 100;\r
63     param.nu = 0.1;\r
64     param.p = 0.3;\r
65 \r
66     param.degree = 1;\r
67     param.coef0 = 0;\r
68     param.shrinking = 1;\r
69     param.probability = 0;\r
70     param.eps = 1e-6;\r
71     param.cache_size = 400;\r
72     param.nr_weight = 0;\r
73     param.weight_label = NULL;\r
74     param.weight = NULL;\r
75     param.kernel_weight = NULL;\r
76     param.kernel_dim = 0;\r
77     param.kernel_norm = 1.;\r
78     param.normalizeKernel = false;\r
79 }\r
80 \r
81 RegressorSVR::~RegressorSVR()\r
82 {\r
83     DEL(node);\r
84 }\r
85 \r
86 struct OptData\r
87 {\r
88     svm_model *svm;\r
89     svm_problem *problem;\r
90 };\r
91 \r
92 double getSVRObjectiveFunction(const svm_model *svm, const double *x, const svm_problem *problem)\r
93 {\r
94     QString gammaString;\r
95     svm_parameter param = svm->param;\r
96     switch(param.kernel_type)\r
97     {\r
98     case LINEAR:\r
99         return 0.;\r
100         break;\r
101     case POLY:\r
102         param.degree = x[0];\r
103         param.gamma = 1. / x[1];\r
104         param.coef0 = x[2];\r
105         break;\r
106     case RBF:\r
107         param.gamma = 1. / x[0];\r
108         break;\r
109     case SIGMOID:\r
110         param.coef0 = x[0];\r
111         break;\r
112     case RBFWEIGH:\r
113     {\r
114         FOR(i, param.kernel_dim)\r
115         {\r
116             param.kernel_weight[i] = x[i];\r
117             gammaString += QString("%1 ").arg(1./x[i]);\r
118         }\r
119     }\r
120         break;\r
121     }\r
122     svm_model *newSVM = svm_train(problem, &param);\r
123     double value = svm_get_dual_objective_function(newSVM);\r
124     qDebug() << "value:" << value << "gamma:" << 1. / param.gamma << "->" << gammaString;\r
125     delete newSVM;\r
126     return value;\r
127 }\r
128 \r
129 double svrObjectiveFunction(unsigned n, const double *x, double *gradient /* NULL if not needed */, void *func_data)\r
130 {\r
131     OptData *data = (OptData*)func_data;\r
132 \r
133     double objective = getSVRObjectiveFunction(data->svm, x, data->problem);\r
134     if(gradient)\r
135     {\r
136         double *dx = new double[n];\r
137         double delta = 1e-2;\r
138         FOR(i, n)\r
139         {\r
140             memcpy(dx, x, n*sizeof(double));\r
141             dx[i] += delta;\r
142             double dError = getSVRObjectiveFunction(data->svm, dx, data->problem);\r
143             gradient[i] = (dError - objective)/delta;\r
144         }\r
145         delete [] dx;\r
146     }\r
147 \r
148     return objective;\r
149 }\r
150 \r
151 void RegressorSVR::Optimize(svm_problem *problem)\r
152 {\r
153     OptData *data = new OptData;\r
154     data->svm = svm;\r
155     data->problem = problem;\r
156     if(svm->param.kernel_type == RBF)\r
157     {\r
158         svm->param.kernel_type = RBFWEIGH;\r
159         svm->param.kernel_weight = new double[dim];\r
160         FOR(d, dim) svm->param.kernel_weight[d] = svm->param.gamma;\r
161         svm->param.gamma = 1.f;\r
162         svm->param.kernel_dim = dim;\r
163     }\r
164 \r
165     int optDim = 1;\r
166     switch(svm->param.kernel_type)\r
167     {\r
168     case POLY:\r
169         optDim = 3;\r
170         break;\r
171     case RBF:\r
172         optDim = 1;\r
173         break;\r
174     case RBFWEIGH:\r
175         optDim = dim;\r
176         break;\r
177     }\r
178 \r
179     nlopt::opt opt(nlopt::LN_AUGLAG, optDim);\r
180     //nlopt::opt opt(nlopt::LN_COBYLA, optDim);\r
181     //nlopt::opt opt(nlopt::LN_NELDERMEAD, optDim);\r
182     //nlopt::opt opt(nlopt::LN_NEWUOA, optDim);\r
183     //nlopt::opt opt(nlopt::LN_PRAXIS, optDim);\r
184     //nlopt::opt opt(nlopt::LN_BOBYQA, optDim);\r
185     //nlopt::opt opt(nlopt::LN_SBPLX, optDim);\r
186 \r
187     opt.set_max_objective(svrObjectiveFunction, (void*)data);\r
188 \r
189     opt.set_maxeval(200);\r
190     vector<double> lowerBounds(optDim, 0.001);\r
191     opt.set_xtol_abs(0.001);\r
192 \r
193     vector<double> x(optDim), xOpt;\r
194 \r
195     vector<double> steps(optDim,0.1);\r
196     switch(svm->param.kernel_type)\r
197     {\r
198     case POLY:\r
199         x[0] = svm->param.degree;\r
200         x[1] = 1. / svm->param.gamma;\r
201         x[2] = svm->param.coef0;\r
202         steps[0] = 1;\r
203         lowerBounds[0] = 1;\r
204         break;\r
205     case RBF:\r
206         x[0] = 1. / svm->param.gamma;\r
207         break;\r
208     case SIGMOID:\r
209         x[0] = svm->param.coef0;\r
210         break;\r
211     case RBFWEIGH:\r
212     {\r
213         FOR(i, svm->param.kernel_dim)\r
214         {\r
215             x[i] = svm->param.kernel_weight[i];\r
216         }\r
217     }\r
218         break;\r
219     }\r
220     opt.set_initial_step(steps);\r
221     opt.set_lower_bounds(lowerBounds);\r
222 \r
223     try\r
224     {\r
225         // do the actual optimization\r
226         xOpt = opt.optimize(x);\r
227         param = svm->param;\r
228         switch(param.kernel_type)\r
229         {\r
230         case POLY:\r
231             param.degree = xOpt[0];\r
232             param.gamma = 1. / xOpt[1];\r
233             param.coef0 = xOpt[2];\r
234             break;\r
235         case RBF:\r
236             param.gamma = 1. / xOpt[0];\r
237             break;\r
238         case SIGMOID:\r
239             param.coef0 = xOpt[0];\r
240             break;\r
241         case RBFWEIGH:\r
242         {\r
243             FOR(i, param.kernel_dim)\r
244             {\r
245                 param.kernel_weight[i] = xOpt[i];\r
246             }\r
247         }\r
248             break;\r
249         }\r
250         delete svm;\r
251         svm = svm_train(problem, &param);\r
252     }\r
253     catch(std::exception e)\r
254     {\r
255         qDebug() << "caught exception while optimizing";\r
256     }\r
257     delete data;\r
258 }\r
259 \r
260 void RegressorSVR::Train(std::vector< fvec > samples, ivec labels)\r
261 {\r
262     svm_problem problem;\r
263     svm_node *x_space;\r
264 \r
265     dim = samples[0].size()-1;\r
266     int oDim = outputDim != -1 && outputDim < dim ? outputDim : dim;\r
267     problem.l = samples.size();\r
268     problem.y = new double[problem.l];\r
269     problem.x = new svm_node *[problem.l];\r
270     x_space = new svm_node[(dim+1)*problem.l];\r
271 \r
272     FOR(i, problem.l)\r
273     {\r
274         FOR(j, dim)\r
275         {\r
276             x_space[(dim+1)*i + j].index = j+1;\r
277             x_space[(dim+1)*i + j].value = samples[i][j];\r
278         }\r
279         x_space[(dim+1)*i + dim].index = -1;\r
280         if(outputDim != -1 && outputDim < dim) x_space[(dim+1)*i + outputDim].value = samples[i][dim];\r
281         problem.x[i] = &x_space[(dim+1)*i];\r
282         problem.y[i] = samples[i][oDim];\r
283     }\r
284 \r
285     DEL(svm);\r
286     DEL(node);\r
287     svm = svm_train(&problem, &param);\r
288     if(bOptimize) Optimize(&problem);\r
289 \r
290     delete [] problem.x;\r
291     delete [] problem.y;\r
292 \r
293     bFixedThreshold = true;\r
294     classThresh = 0.5f;\r
295 }\r
296 \r
297 fvec RegressorSVR::Test( const fvec &sample )\r
298 {\r
299     int dim = sample.size()-1;\r
300     float estimate;\r
301     if(!node)\r
302     {\r
303         node = new svm_node[dim+1];\r
304         node[dim].index = -1;\r
305     }\r
306     FOR(i, dim)\r
307     {\r
308         node[i].index = i+1;\r
309         node[i].value = sample[i];\r
310     }\r
311     if(outputDim != -1 && outputDim < dim) node[outputDim].value = sample[dim];\r
312     estimate = (float)svm_predict(svm, node);\r
313     fvec res;\r
314     res.push_back(estimate);\r
315     res.push_back(1);\r
316     return res;\r
317 }\r
318 \r
319 fVec RegressorSVR::Test( const fVec &sample )\r
320 {\r
321     int dim = 1;\r
322     float estimate;\r
323     if(!node) node = new svm_node[dim+1];\r
324     FOR(i, dim)\r
325     {\r
326         node[i].index = i+1;\r
327         node[i].value = sample._[i];\r
328     }\r
329     node[dim].index = -1;\r
330     estimate = (float)svm_predict(svm, node);\r
331     return fVec(estimate,1);\r
332 }\r
333 \r
334 void RegressorSVR::SetParams(int svmType, float svmC, float svmP, u32 kernelType, float kernelParam)\r
335 {\r
336     // default values\r
337     param.svm_type = svmType;\r
338     param.C = svmC;\r
339     param.nu = svmC;\r
340     param.eps = 0.01;\r
341     param.p = svmP;\r
342 \r
343     param.coef0 = 0;\r
344     param.gamma = 1;\r
345 \r
346     switch(kernelType)\r
347     {\r
348     case 0:\r
349         param.kernel_type = LINEAR;\r
350         param.degree = 1;\r
351         break;\r
352     case 1:\r
353         param.kernel_type = POLY;\r
354         param.degree = (u32)kernelParam;\r
355         break;\r
356     case 2:\r
357         param.kernel_type = RBF;\r
358         param.gamma = kernelParam;\r
359         break;\r
360     case 3:\r
361         param.kernel_type = SIGMOID;\r
362         param.gamma = kernelParam;\r
363         break;\r
364     }\r
365 }\r