Aluode/PerceptionLabPortable
0
1import numpy as np
2from scipy import sparse
3from ..utils._cython_blas cimport _dot
4from ..utils._typedefs cimport float64_t, int32_t, intp_t
5
6cdef extern from *:
7 ctypedef char* const_char_p "const char*"
8
9################################################################################
10# Includes
11
12cdef extern from "_svm_cython_blas_helpers.h":
13 ctypedef double (*dot_func)(int, const double*, int, const double*, int)
14 cdef struct BlasFunctions:
15 dot_func dot
16
17cdef extern from "svm.h":
18 cdef struct svm_csr_node
19 cdef struct svm_csr_model
20 cdef struct svm_parameter
21 cdef struct svm_csr_problem
22 char *svm_csr_check_parameter(svm_csr_problem *, svm_parameter *)
23 svm_csr_model *svm_csr_train(svm_csr_problem *, svm_parameter *, int *, BlasFunctions *) nogil
24 void svm_csr_free_and_destroy_model(svm_csr_model** model_ptr_ptr)
25
26cdef extern from "libsvm_sparse_helper.c":
27 # this file contains methods for accessing libsvm 'hidden' fields
28 svm_csr_problem * csr_set_problem (
29 char *, intp_t *, char *, intp_t *, char *, char *, char *, int)
30 svm_csr_model *csr_set_model(svm_parameter *param, int nr_class,
31 char *SV_data, intp_t *SV_indices_dims,
32 char *SV_indices, intp_t *SV_intptr_dims,
33 char *SV_intptr,
34 char *sv_coef, char *rho, char *nSV,
35 char *probA, char *probB)
36 svm_parameter *set_parameter (int , int , int , double, double ,
37 double , double , double , double,
38 double, int, int, int, char *, char *, int,
39 int)
40 void copy_sv_coef (char *, svm_csr_model *)
41 void copy_n_iter (char *, svm_csr_model *)
42 void copy_support (char *, svm_csr_model *)
43 void copy_intercept (char *, svm_csr_model *, intp_t *)
44 int copy_predict (char *, svm_csr_model *, intp_t *, char *, BlasFunctions *)
45 int csr_copy_predict_values (intp_t *data_size, char *data, intp_t *index_size,
46 char *index, intp_t *intptr_size, char *size,
47 svm_csr_model *model, char *dec_values, int nr_class, BlasFunctions *)
48 int csr_copy_predict (intp_t *data_size, char *data, intp_t *index_size,
49 char *index, intp_t *intptr_size, char *size,
50 svm_csr_model *model, char *dec_values, BlasFunctions *) nogil
51 int csr_copy_predict_proba (intp_t *data_size, char *data, intp_t *index_size,
52 char *index, intp_t *intptr_size, char *size,
53 svm_csr_model *model, char *dec_values, BlasFunctions *) nogil
54
55 int copy_predict_values(char *, svm_csr_model *, intp_t *, char *, int, BlasFunctions *)
56 int csr_copy_SV (char *values, intp_t *n_indices,
57 char *indices, intp_t *n_indptr, char *indptr,
58 svm_csr_model *model, int n_features)
59 intp_t get_nonzero_SV (svm_csr_model *)
60 void copy_nSV (char *, svm_csr_model *)
61 void copy_probA (char *, svm_csr_model *, intp_t *)
62 void copy_probB (char *, svm_csr_model *, intp_t *)
63 intp_t get_l (svm_csr_model *)
64 intp_t get_nr (svm_csr_model *)
65 int free_problem (svm_csr_problem *)
66 int free_model (svm_csr_model *)
67 int free_param (svm_parameter *)
68 int free_model_SV(svm_csr_model *model)
69 void set_verbosity(int)
70
71
72def libsvm_sparse_train (int n_features,
73 const float64_t[::1] values,
74 const int32_t[::1] indices,
75 const int32_t[::1] indptr,
76 const float64_t[::1] Y,
77 int svm_type, int kernel_type, int degree, double gamma,
78 double coef0, double eps, double C,
79 const float64_t[::1] class_weight,
80 const float64_t[::1] sample_weight,
81 double nu, double cache_size, double p, int
82 shrinking, int probability, int max_iter,
83 int random_seed):
84 """
85 Wrap svm_train from libsvm using a scipy.sparse.csr matrix
86
87 Work in progress.
88
89 Parameters
90 ----------
91 n_features : number of features.
92 XXX: can we retrieve this from any other parameter ?
93
94 X : array-like, dtype=float, size=[N, D]
95
96 Y : array, dtype=float, size=[N]
97 target vector
98
99 ...
100
101 Notes
102 -------------------
103 See sklearn.svm.predict for a complete list of parameters.
104
105 """
106
107 cdef svm_parameter *param
108 cdef svm_csr_problem *problem
109 cdef svm_csr_model *model
110 cdef const_char_p error_msg
111
112 if len(sample_weight) == 0:
113 sample_weight = np.ones(Y.shape[0], dtype=np.float64)
114 else:
115 assert sample_weight.shape[0] == indptr.shape[0] - 1, \
116 "sample_weight and X have incompatible shapes: " + \
117 "sample_weight has %s samples while X has %s" % \
118 (sample_weight.shape[0], indptr.shape[0] - 1)
119
120 # we should never end up here with a precomputed kernel matrix,
121 # as this is always dense.
122 assert(kernel_type != 4)
123
124 # set libsvm problem
125 problem = csr_set_problem(
126 <char *> &values[0],
127 <intp_t *> indices.shape,
128 <char *> &indices[0],
129 <intp_t *> indptr.shape,
130 <char *> &indptr[0],
131 <char *> &Y[0],
132 <char *> &sample_weight[0],
133 kernel_type,
134 )
135
136 cdef int32_t[::1] \
137 class_weight_label = np.arange(class_weight.shape[0], dtype=np.int32)
138
139 # set parameters
140 param = set_parameter(
141 svm_type,
142 kernel_type,
143 degree,
144 gamma,
145 coef0,
146 nu,
147 cache_size,
148 C,
149 eps,
150 p,
151 shrinking,
152 probability,
153 <int> class_weight.shape[0],
154 <char *> &class_weight_label[0] if class_weight_label.size > 0 else NULL,
155 <char *> &class_weight[0] if class_weight.size > 0 else NULL, max_iter,
156 random_seed,
157 )
158
159 # check parameters
160 if (param == NULL or problem == NULL):
161 raise MemoryError("Seems we've run out of memory")
162 error_msg = svm_csr_check_parameter(problem, param)
163 if error_msg:
164 free_problem(problem)
165 free_param(param)
166 raise ValueError(error_msg)
167 cdef BlasFunctions blas_functions
168 blas_functions.dot = _dot[double]
169 # call svm_train, this does the real work
170 cdef int fit_status = 0
171 with nogil:
172 model = svm_csr_train(problem, param, &fit_status, &blas_functions)
173
174 cdef intp_t SV_len = get_l(model)
175 cdef intp_t n_class = get_nr(model)
176
177 cdef int[::1] n_iter
178 n_iter = np.empty(max(1, n_class * (n_class - 1) // 2), dtype=np.intc)
179 copy_n_iter(<char *> &n_iter[0], model)
180
181 # copy model.sv_coef
182 # we create a new array instead of resizing, otherwise
183 # it would not erase previous information
184 cdef float64_t[::1] sv_coef_data
185 sv_coef_data = np.empty((n_class-1)*SV_len, dtype=np.float64)
186 copy_sv_coef (<char *> &sv_coef_data[0] if sv_coef_data.size > 0 else NULL, model)
187
188 cdef int32_t[::1] support
189 support = np.empty(SV_len, dtype=np.int32)
190 copy_support(<char *> &support[0] if support.size > 0 else NULL, model)
191
192 # copy model.rho into the intercept
193 # the intercept is just model.rho but with sign changed
194 cdef float64_t[::1]intercept
195 intercept = np.empty(n_class*(n_class-1)//2, dtype=np.float64)
196 copy_intercept (<char *> &intercept[0], model, <intp_t *> intercept.shape)
197
198 # copy model.SV
199 # we erase any previous information in SV
200 # TODO: custom kernel
201 cdef intp_t nonzero_SV
202 nonzero_SV = get_nonzero_SV (model)
203
204 cdef float64_t[::1] SV_data
205 cdef int32_t[::1] SV_indices, SV_indptr
206 SV_data = np.empty(nonzero_SV, dtype=np.float64)
207 SV_indices = np.empty(nonzero_SV, dtype=np.int32)
208 SV_indptr = np.empty(<intp_t>SV_len + 1, dtype=np.int32)
209 csr_copy_SV(
210 <char *> &SV_data[0] if SV_data.size > 0 else NULL,
211 <intp_t *> SV_indices.shape,
212 <char *> &SV_indices[0] if SV_indices.size > 0 else NULL,
213 <intp_t *> SV_indptr.shape,
214 <char *> &SV_indptr[0] if SV_indptr.size > 0 else NULL,
215 model,
216 n_features,
217 )
218 support_vectors_ = sparse.csr_matrix(
219 (SV_data, SV_indices, SV_indptr), (SV_len, n_features)
220 )
221
222 # copy model.nSV
223 # TODO: do only in classification
224 cdef int32_t[::1]n_class_SV
225 n_class_SV = np.empty(n_class, dtype=np.int32)
226 copy_nSV(<char *> &n_class_SV[0], model)
227
228 # # copy probabilities
229 cdef float64_t[::1] probA, probB
230 if probability != 0:
231 if svm_type < 2: # SVC and NuSVC
232 probA = np.empty(n_class*(n_class-1)//2, dtype=np.float64)
233 probB = np.empty(n_class*(n_class-1)//2, dtype=np.float64)
234 copy_probB(<char *> &probB[0], model, <intp_t *> probB.shape)
235 else:
236 probA = np.empty(1, dtype=np.float64)
237 probB = np.empty(0, dtype=np.float64)
238 copy_probA(<char *> &probA[0], model, <intp_t *> probA.shape)
239 else:
240 probA = np.empty(0, dtype=np.float64)
241 probB = np.empty(0, dtype=np.float64)
242
243 svm_csr_free_and_destroy_model (&model)
244 free_problem(problem)
245 free_param(param)
246
247 return (
248 support.base,
249 support_vectors_,
250 sv_coef_data.base,
251 intercept.base,
252 n_class_SV.base,
253 probA.base,
254 probB.base,
255 fit_status,
256 n_iter.base,
257 )
258
259
260def libsvm_sparse_predict (const float64_t[::1] T_data,
261 const int32_t[::1] T_indices,
262 const int32_t[::1] T_indptr,
263 const float64_t[::1] SV_data,
264 const int32_t[::1] SV_indices,
265 const int32_t[::1] SV_indptr,
266 const float64_t[::1] sv_coef,
267 const float64_t[::1]
268 intercept, int svm_type, int kernel_type, int
269 degree, double gamma, double coef0, double
270 eps, double C,
271 const float64_t[:] class_weight,
272 double nu, double p, int
273 shrinking, int probability,
274 const int32_t[::1] nSV,
275 const float64_t[::1] probA,
276 const float64_t[::1] probB):
277 """
278 Predict values T given a model.
279
280 For speed, all real work is done at the C level in function
281 copy_predict (libsvm_helper.c).
282
283 We have to reconstruct model and parameters to make sure we stay
284 in sync with the python object.
285
286 See sklearn.svm.predict for a complete list of parameters.
287
288 Parameters
289 ----------
290 X : array-like, dtype=float
291 Y : array
292 target vector
293
294 Returns
295 -------
296 dec_values : array
297 predicted values.
298 """
299 cdef float64_t[::1] dec_values
300 cdef svm_parameter *param
301 cdef svm_csr_model *model
302 cdef int32_t[::1] \
303 class_weight_label = np.arange(class_weight.shape[0], dtype=np.int32)
304 cdef int rv
305 param = set_parameter(
306 svm_type,
307 kernel_type,
308 degree,
309 gamma,
310 coef0,
311 nu,
312 100.0, # cache size has no effect on predict
313 C,
314 eps,
315 p,
316 shrinking,
317 probability,
318 <int> class_weight.shape[0],
319 <char *> &class_weight_label[0] if class_weight_label.size > 0 else NULL,
320 <char *> &class_weight[0] if class_weight.size > 0 else NULL,
321 -1,
322 -1, # random seed has no effect on predict either
323 )
324
325 model = csr_set_model(
326 param, <int> nSV.shape[0],
327 <char *> &SV_data[0] if SV_data.size > 0 else NULL,
328 <intp_t *>SV_indices.shape,
329 <char *> &SV_indices[0] if SV_indices.size > 0 else NULL,
330 <intp_t *> SV_indptr.shape,
331 <char *> &SV_indptr[0] if SV_indptr.size > 0 else NULL,
332 <char *> &sv_coef[0] if sv_coef.size > 0 else NULL,
333 <char *> &intercept[0],
334 <char *> &nSV[0],
335 <char *> &probA[0] if probA.size > 0 else NULL,
336 <char *> &probB[0] if probB.size > 0 else NULL,
337 )
338 # TODO: use check_model
339 dec_values = np.empty(T_indptr.shape[0]-1)
340 cdef BlasFunctions blas_functions
341 blas_functions.dot = _dot[double]
342 with nogil:
343 rv = csr_copy_predict(
344 <intp_t *> T_data.shape,
345 <char *> &T_data[0],
346 <intp_t *> T_indices.shape,
347 <char *> &T_indices[0],
348 <intp_t *> T_indptr.shape,
349 <char *> &T_indptr[0],
350 model,
351 <char *> &dec_values[0],
352 &blas_functions,
353 )
354 if rv < 0:
355 raise MemoryError("We've run out of memory")
356 # free model and param
357 free_model_SV(model)
358 free_model(model)
359 free_param(param)
360 return dec_values.base
361
362
363def libsvm_sparse_predict_proba(
364 const float64_t[::1] T_data,
365 const int32_t[::1] T_indices,
366 const int32_t[::1] T_indptr,
367 const float64_t[::1] SV_data,
368 const int32_t[::1] SV_indices,
369 const int32_t[::1] SV_indptr,
370 const float64_t[::1] sv_coef,
371 const float64_t[::1]
372 intercept, int svm_type, int kernel_type, int
373 degree, double gamma, double coef0, double
374 eps, double C,
375 const float64_t[:] class_weight,
376 double nu, double p, int shrinking, int probability,
377 const int32_t[::1] nSV,
378 const float64_t[::1] probA,
379 const float64_t[::1] probB,
380):
381 """
382 Predict values T given a model.
383 """
384 cdef float64_t[:, ::1] dec_values
385 cdef svm_parameter *param
386 cdef svm_csr_model *model
387 cdef int32_t[::1] \
388 class_weight_label = np.arange(class_weight.shape[0], dtype=np.int32)
389 param = set_parameter(
390 svm_type,
391 kernel_type,
392 degree,
393 gamma,
394 coef0,
395 nu,
396 100.0, # cache size has no effect on predict
397 C,
398 eps,
399 p,
400 shrinking,
401 probability,
402 <int> class_weight.shape[0],
403 <char *> &class_weight_label[0] if class_weight_label.size > 0 else NULL,
404 <char *> &class_weight[0] if class_weight.size > 0 else NULL,
405 -1,
406 -1, # random seed has no effect on predict either
407 )
408
409 model = csr_set_model(
410 param,
411 <int> nSV.shape[0],
412 <char *> &SV_data[0] if SV_data.size > 0 else NULL,
413 <intp_t *> SV_indices.shape,
414 <char *> &SV_indices[0] if SV_indices.size > 0 else NULL,
415 <intp_t *> SV_indptr.shape,
416 <char *> &SV_indptr[0] if SV_indptr.size > 0 else NULL,
417 <char *> &sv_coef[0] if sv_coef.size > 0 else NULL,
418 <char *> &intercept[0],
419 <char *> &nSV[0],
420 <char *> &probA[0] if probA.size > 0 else NULL,
421 <char *> &probB[0] if probB.size > 0 else NULL,
422 )
423 # TODO: use check_model
424 cdef intp_t n_class = get_nr(model)
425 cdef int rv
426 dec_values = np.empty((T_indptr.shape[0]-1, n_class), dtype=np.float64)
427 cdef BlasFunctions blas_functions
428 blas_functions.dot = _dot[double]
429 with nogil:
430 rv = csr_copy_predict_proba(
431 <intp_t *> T_data.shape,
432 <char *> &T_data[0],
433 <intp_t *> T_indices.shape,
434 <char *> &T_indices[0],
435 <intp_t *> T_indptr.shape,
436 <char *> &T_indptr[0],
437 model,
438 <char *> &dec_values[0, 0],
439 &blas_functions,
440 )
441 if rv < 0:
442 raise MemoryError("We've run out of memory")
443 # free model and param
444 free_model_SV(model)
445 free_model(model)
446 free_param(param)
447 return dec_values.base
448
449
450def libsvm_sparse_decision_function(
451 const float64_t[::1] T_data,
452 const int32_t[::1] T_indices,
453 const int32_t[::1] T_indptr,
454 const float64_t[::1] SV_data,
455 const int32_t[::1] SV_indices,
456 const int32_t[::1] SV_indptr,
457 const float64_t[::1] sv_coef,
458 const float64_t[::1]
459 intercept, int svm_type, int kernel_type, int
460 degree, double gamma, double coef0, double
461 eps, double C,
462 const float64_t[:] class_weight,
463 double nu, double p, int shrinking, int probability,
464 const int32_t[::1] nSV,
465 const float64_t[::1] probA,
466 const float64_t[::1] probB,
467):
468 """
469 Predict margin (libsvm name for this is predict_values)
470
471 We have to reconstruct model and parameters to make sure we stay
472 in sync with the python object.
473 """
474 cdef float64_t[:, ::1] dec_values
475 cdef svm_parameter *param
476 cdef intp_t n_class
477
478 cdef svm_csr_model *model
479 cdef int32_t[::1] \
480 class_weight_label = np.arange(class_weight.shape[0], dtype=np.int32)
481 param = set_parameter(
482 svm_type,
483 kernel_type,
484 degree,
485 gamma,
486 coef0,
487 nu,
488 100.0, # cache size has no effect on predict
489 C,
490 eps,
491 p,
492 shrinking,
493 probability,
494 <int> class_weight.shape[0],
495 <char *> &class_weight_label[0] if class_weight_label.size > 0 else NULL,
496 <char *> &class_weight[0] if class_weight.size > 0 else NULL,
497 -1,
498 -1,
499 )
500
501 model = csr_set_model(
502 param,
503 <int> nSV.shape[0],
504 <char *> &SV_data[0] if SV_data.size > 0 else NULL,
505 <intp_t *> SV_indices.shape,
506 <char *> &SV_indices[0] if SV_indices.size > 0 else NULL,
507 <intp_t *> SV_indptr.shape,
508 <char *> &SV_indptr[0] if SV_indptr.size > 0 else NULL,
509 <char *> &sv_coef[0] if sv_coef.size > 0 else NULL,
510 <char *> &intercept[0],
511 <char *> &nSV[0],
512 <char *> &probA[0] if probA.size > 0 else NULL,
513 <char *> &probB[0] if probB.size > 0 else NULL,
514 )
515
516 if svm_type > 1:
517 n_class = 1
518 else:
519 n_class = get_nr(model)
520 n_class = n_class * (n_class - 1) // 2
521
522 dec_values = np.empty((T_indptr.shape[0] - 1, n_class), dtype=np.float64)
523 cdef BlasFunctions blas_functions
524 blas_functions.dot = _dot[double]
525 if csr_copy_predict_values(
526 <intp_t *> T_data.shape,
527 <char *> &T_data[0],
528 <intp_t *> T_indices.shape,
529 <char *> &T_indices[0],
530 <intp_t *> T_indptr.shape,
531 <char *> &T_indptr[0],
532 model,
533 <char *> &dec_values[0, 0],
534 n_class,
535 &blas_functions,
536 ) < 0:
537 raise MemoryError("We've run out of memory")
538 # free model and param
539 free_model_SV(model)
540 free_model(model)
541 free_param(param)
542
543 return dec_values.base
544
545
546def set_verbosity_wrap(int verbosity):
547 """
548 Control verbosity of libsvm library
549 """
550 set_verbosity(verbosity)
551