main.cpp 28.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
/**
 * @author Andre Anjos <andre.anjos@idiap.ch>
 * @date Tue  3 Dec 14:13:22 2013 CET
 *
 * @brief Bindings to bob::math
 */

#ifdef NO_IMPORT_ARRAY
#undef NO_IMPORT_ARRAY
#endif
#include <xbob.blitz/capi.h>
André Anjos's avatar
André Anjos committed
12
#include <xbob.blitz/cleanup.h>
13

14
15
#include <xbob.extension/documentation.h>

16
#include "histogram.h"
André Anjos's avatar
André Anjos committed
17
#include "linsolve.h"
André Anjos's avatar
André Anjos committed
18
#include "pavx.h"
André Anjos's avatar
André Anjos committed
19
#include "norminv.h"
20
#include "scatter.h"
21
#include "lp_interior_point.h"
22

23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
static xbob::extension::FunctionDoc s_histogram_intersection = xbob::extension::FunctionDoc(
    "histogram_intersection",
    "Computes the histogram intersection between the given histograms, which might be of singular dimension only.",
    "The histogram intersection is computed as follows:\n"
    ".. math:: sim(h_1,h_2) = \\sum_i \\min \\{h_{1i}, h_{2i}\\}\n"
    "The histogram intersection defines a similarity measure, so higher values are better. "
    "You can use this method in two different formats. "
    "The first interface accepts non-sparse histograms. "
    "The second interface accepts sparse histograms represented by indexes and values.\n"
    ".. note:: Histograms are given as two matrices, one with the indexes and one with the data. All data points that for which no index exists are considered to be zero.\n"
    ".. note:: In general, histogram intersection with sparse histograms needs more time to be computed."
  )
  .add_prototype("h1, h2", "sim")
  .add_prototype("index_1, value_1, index_2, value_2", "sim")
  .add_parameter("h1, h2", "array_like (1D)", "Histograms to compute the histogram intersection for")
  .add_parameter("index_1, index_2", "array_like (int, 1D)", "Indices of the sparse histograms value_1 and value_2")
  .add_parameter("value_1, value_2", "array_like (1D)", "Sparse histograms to compute the histogram intersection for")
  .add_return("sim", "float", "The histogram intersection value for the given histograms.")
;

static xbob::extension::FunctionDoc s_chi_square = xbob::extension::FunctionDoc(
    "chi_square",
    "Computes the chi square distance between the given histograms, which might be of singular dimension only.",
    "The chi square distance is computed as follows:\n"
    ".. math:: dist(h_1,h_2) = \\sum_i \\frac{(h_{1i} - h_{2i})^2}{h_{1i} + h_{2i}}\n"
    "Chi square defines a distance metric, so lower values are better. "
    "You can use this method in two different formats. "
    "The first interface accepts non-sparse histograms. "
    "The second interface accepts sparse histograms represented by indexes and values.\n"
    ".. note:: Histograms are given as two matrices, one with the indexes and one with the data. All data points that for which no index exists are considered to be zero.\n"
    ".. note:: In general, histogram intersection with sparse histograms needs more time to be computed."
  )
  .add_prototype("h1, h2", "dist")
  .add_prototype("index_1, value_1, index_2, value_2", "dist")
  .add_parameter("h1, h2", "array_like (1D)", "Histograms to compute the chi square distance for")
  .add_parameter("index_1, index_2", "array_like (int, 1D)", "Indices of the sparse histograms value_1 and value_2")
  .add_parameter("value_1, value_2", "array_like (1D)", "Sparse histograms to compute the chi square distance for")
  .add_return("dist", "float", "The chi square distance value for the given histograms.")
;

static xbob::extension::FunctionDoc s_kullback_leibler = xbob::extension::FunctionDoc(
    "kullback_leibler",
    "Computes the Kullback-Leibler histogram divergence between the given histograms, which might be of singular dimension only.",
    "The chi square distance is inspired by `link <http://www.informatik.uni-freiburg.de/~tipaldi/FLIRTLib/HistogramDistances_8hpp_source.html>`_ and computed as follows:\n"
    ".. math:: dist(h_1,h_2) = \\sum_i (h_{1i} - h_{2i}) * \\log (h_{1i} / h_{2i})\n"
    "The Kullback-Leibler divergence defines a distance metric, so lower values are better. "
    "You can use this method in two different formats. "
    "The first interface accepts non-sparse histograms. "
    "The second interface accepts sparse histograms represented by indexes and values.\n"
    ".. note:: Histograms are given as two matrices, one with the indexes and one with the data. All data points that for which no index exists are considered to be zero.\n"
    ".. note:: In general, histogram intersection with sparse histograms needs more time to be computed."
  )
  .add_prototype("h1, h2", "dist")
  .add_prototype("index_1, value_1, index_2, value_2", "dist")
  .add_parameter("h1, h2", "array_like (1D)", "Histograms to compute the Kullback-Leibler divergence for")
  .add_parameter("index_1, index_2", "array_like (int, 1D)", "Indices of the sparse histograms value_1 and value_2")
  .add_parameter("value_1, value_2", "array_like (1D)", "Sparse histograms to compute the Kullback-Leibler divergence for")
  .add_return("dist", "float", "The Kullback-Leibler divergence value for the given histograms.")
;
82

83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
static xbob::extension::FunctionDoc s_linsolve = xbob::extension::FunctionDoc(
  "linsolve",
  "Solves the linear system :math:`Ax=b` and returns the result in :math:`x`.",
  "This method uses LAPACK's ``dgesv`` generic solver. "
  "You can use this method in two different formats. "
  "The first interface accepts the matrices :math:`A` and :math:`b` returning :math:`x`. "
  "The second one accepts a pre-allocated :math:`x` vector and sets it with the linear system solution."
  )
  .add_prototype("A, b", "x")
  .add_prototype("A, b, x")
  .add_parameter("A", "array_like (2D)", "The matrix :math:`A` of the linear system")
  .add_parameter("b", "array_like (1D)", "The vector :math:`b` of the linear system")
  .add_parameter("x", "array_like (1D)", "The result vector :math:`x`, as parameter")
  .add_return("x", "array_like (1D)", "The result vector :math:`x`, as return value")
;
98

99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
static xbob::extension::FunctionDoc s_linsolve_nocheck = xbob::extension::FunctionDoc(
  "linsolve_",
  "Solves the linear system :math:`Ax=b` and returns the result in :math:`x`.",
  ".. warning:: This variant does not perform any checks on the input matrices and is faster then :func:`linsolve`. "
  "Use it when you are sure your input matrices sizes match.\n"
  "This method uses LAPACK's ``dgesv`` generic solver. "
  "You can use this method in two different formats. "
  "The first interface accepts the matrices :math:`A` and :math:`b` returning :math:`x`. "
  "The second one accepts a pre-allocated :math:`x` vector and sets it with the linear system solution.\n"
  )
  .add_prototype("A, b", "x")
  .add_prototype("A, b, x")
  .add_parameter("A", "array_like (2D)", "The matrix :math:`A` of the linear system")
  .add_parameter("b", "array_like (1D)", "The vector :math:`b` of the linear system")
  .add_parameter("x", "array_like (1D)", "The result vector :math:`x`, as parameter")
  .add_return("x", "array_like (1D)", "The result vector :math:`x`, as return value")
;
116

117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
static xbob::extension::FunctionDoc s_linsolve_sympos = xbob::extension::FunctionDoc(
  "linsolve_sympos",
  "Solves the linear system :math:`Ax=b` and returns the result in :math:`x` for symmetric :math:`A` matrix.",
  "This method uses LAPACK's ``dposv`` solver, assuming :math:`A` is a symmetric positive definite matrix. "
  "You can use this method in two different formats. "
  "The first interface accepts the matrices :math:`A` and :math:`b` returning :math:`x`. "
  "The second one accepts a pre-allocated :math:`x` vector and sets it with the linear system solution."
  )
  .add_prototype("A, b", "x")
  .add_prototype("A, b, x")
  .add_parameter("A", "array_like (2D)", "The matrix :math:`A` of the linear system")
  .add_parameter("b", "array_like (1D)", "The vector :math:`b` of the linear system")
  .add_parameter("x", "array_like (1D)", "The result vector :math:`x`, as parameter")
  .add_return("x", "array_like (1D)", "The result vector :math:`x`, as return value")
;
André Anjos's avatar
André Anjos committed
132

133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
static xbob::extension::FunctionDoc s_linsolve_sympos_nocheck = xbob::extension::FunctionDoc(
  "linsolve_sympos_",
  ".. warning:: This variant does not perform any checks on the input matrices and is faster then :func:`linsolve_sympos`. "
  "Use it when you are sure your input matrices sizes match.\n"
  "Solves the linear system :math:`Ax=b` and returns the result in :math:`x` for symmetric :math:`A` matrix.",
  "This method uses LAPACK's ``dposv`` solver, assuming :math:`A` is a symmetric positive definite matrix. "
  "You can use this method in two different formats. "
  "The first interface accepts the matrices :math:`A` and :math:`b` returning :math:`x`. "
  "The second one accepts a pre-allocated :math:`x` vector and sets it with the linear system solution."
  )
  .add_prototype("A, b", "x")
  .add_prototype("A, b, x")
  .add_parameter("A", "array_like (2D)", "The matrix :math:`A` of the linear system")
  .add_parameter("b", "array_like (1D)", "The vector :math:`b` of the linear system")
  .add_parameter("x", "array_like (1D)", "The result vector :math:`x`, as parameter")
  .add_return("x", "array_like (1D)", "The result vector :math:`x`, as return value")
;
André Anjos's avatar
André Anjos committed
150

151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
static xbob::extension::FunctionDoc s_linsolve_cg_sympos = xbob::extension::FunctionDoc(
  "linsolve_cg_sympos",
  "Solves the linear system :math:`Ax=b` using conjugate gradients and returns the result in :math:`x` for symmetric :math:`A` matrix.",
  "This method uses the conjugate gradient solver, assuming :math:`A` is a symmetric positive definite matrix. "
  "You can use this method in two different formats. "
  "The first interface accepts the matrices :math:`A` and :math:`b` returning :math:`x`. "
  "The second one accepts a pre-allocated :math:`x` vector and sets it with the linear system solution."
  )
  .add_prototype("A, b", "x")
  .add_prototype("A, b, x")
  .add_parameter("A", "array_like (2D)", "The matrix :math:`A` of the linear system")
  .add_parameter("b", "array_like (1D)", "The vector :math:`b` of the linear system")
  .add_parameter("x", "array_like (1D)", "The result vector :math:`x`, as parameter")
  .add_return("x", "array_like (1D)", "The result vector :math:`x`, as return value")
;
André Anjos's avatar
André Anjos committed
166

167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
static xbob::extension::FunctionDoc s_linsolve_cg_sympos_nocheck = xbob::extension::FunctionDoc(
  "linsolve_cg_sympos_",
  "Solves the linear system :math:`Ax=b` using conjugate gradients and returns the result in :math:`x` for symmetric :math:`A` matrix.",
  ".. warning:: This variant does not perform any checks on the input matrices and is faster then :func:`linsolve_cg_sympos`. "
  "Use it when you are sure your input matrices sizes match.\n"
  "This method uses the conjugate gradient solver, assuming :math:`A` is a symmetric positive definite matrix. "
  "You can use this method in two different formats. "
  "The first interface accepts the matrices :math:`A` and :math:`b` returning :math:`x`. "
  "The second one accepts a pre-allocated :math:`x` vector and sets it with the linear system solution."
  )
  .add_prototype("A, b", "x")
  .add_prototype("A, b, x")
  .add_parameter("A", "array_like (2D)", "The matrix :math:`A` of the linear system")
  .add_parameter("b", "array_like (1D)", "The vector :math:`b` of the linear system")
  .add_parameter("x", "array_like (1D)", "The result vector :math:`x`, as parameter")
  .add_return("x", "array_like (1D)", "The result vector :math:`x`, as return value")
;
André Anjos's avatar
André Anjos committed
184

185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
static xbob::extension::FunctionDoc s_pavx = xbob::extension::FunctionDoc(
  "pavx",
  "Applies the Pool-Adjacent-Violators Algorithm",
  "Applies the Pool-Adjacent-Violators Algorithm to ``input``. "
  "This is a simplified C++ port of the isotonic regression code made available at the `University of Bern website <http://www.imsv.unibe.ch/content/staff/personalhomepages/duembgen/software/isotonicregression/index_eng.html>`_.\n"
  "You can use this method in two different formats. "
  "The first interface accepts the ``input`` and ``output``. "
  "The second one accepts the input array ``input`` and allocates a new ``output`` array, which is returned. "
  )
  .add_prototype("input, output")
  .add_prototype("input", "output")
  .add_parameter("input", "array_like (float, 1D)", "The input matrix for the PAV algorithm.")
  .add_parameter("output", "array_like (float, 1D)", "The output matrix, must be of the same size as ``input``")
  .add_return("output", "array_like (float, 1D)", "The output matrix; will be created in the same size as ``input``")
;
André Anjos's avatar
André Anjos committed
200

201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
static xbob::extension::FunctionDoc s_pavx_nocheck = xbob::extension::FunctionDoc(
  "pavx_",
  "Applies the Pool-Adjacent-Violators Algorithm",
  ".. warning:: This variant does not perform any checks on the input matrices and is faster then :func:`pavx`. "
  "Use it when you are sure your input matrices sizes match.\n"
  "Applies the Pool-Adjacent-Violators Algorithm to ``input``. "
  "This is a simplified C++ port of the isotonic regression code made available at the `University of Bern website <http://www.imsv.unibe.ch/content/staff/personalhomepages/duembgen/software/isotonicregression/index_eng.html>`_.\n"
  "You can use this method in two different formats. "
  "The first interface accepts the ``input`` and ``output``. "
  "The second one accepts the input array ``input`` and allocates a new ``output`` array, which is returned. "
  )
  .add_prototype("input, output")
  .add_prototype("input", "output")
  .add_parameter("input", "array_like (float, 1D)", "The input matrix for the PAV algorithm.")
  .add_parameter("output", "array_like (float, 1D)", "The output matrix, must be of the same size as ``input``")
  .add_return("output", "array_like (float, 1D)", "The output matrix; will be created in the same size as ``input``")
;
André Anjos's avatar
André Anjos committed
218

219
220
221
222
223
224
225
226
227
228
229
static xbob::extension::FunctionDoc s_pavx_width = xbob::extension::FunctionDoc(
  "pavxWidth",
  "Applies the Pool-Adjacent-Violators Algorithm and returns the width.",
  "Applies the Pool-Adjacent-Violators Algorithm to ``input``. "
  "This is a simplified C++ port of the isotonic regression code made available at the `University of Bern website <http://www.imsv.unibe.ch/content/staff/personalhomepages/duembgen/software/isotonicregression/index_eng.html>`_."
  )
  .add_prototype("input, output", "width")
  .add_parameter("input", "array_like (float, 1D)", "The input matrix for the PAV algorithm.")
  .add_parameter("output", "array_like (float, 1D)", "The output matrix, must be of the same size as ``input``")
  .add_return("width", "array_like (uint64, 1D)", "The width matrix will be created in the same size as ``input``\n.. todo:: Explain, what width means in this case")
;
André Anjos's avatar
André Anjos committed
230

231
232
233
234
235
236
237
238
239
240
241
242
static xbob::extension::FunctionDoc s_pavx_width_height = xbob::extension::FunctionDoc(
  "pavxWidthHeight",
  "Applies the Pool-Adjacent-Violators Algorithm and returns the width and the height.",
  "Applies the Pool-Adjacent-Violators Algorithm to ``input``. "
  "This is a simplified C++ port of the isotonic regression code made available at the `University of Bern website <http://www.imsv.unibe.ch/content/staff/personalhomepages/duembgen/software/isotonicregression/index_eng.html>`_."
  )
  .add_prototype("input, output", "width, height")
  .add_parameter("input", "array_like (float, 1D)", "The input matrix for the PAV algorithm.")
  .add_parameter("output", "array_like (float, 1D)", "The output matrix, must be of the same size as ``input``")
  .add_return("width", "array_like (uint64, 1D)", "The width matrix will be created in the same size as ``input``\n.. todo:: Explain, what width means in this case")
  .add_return("height", "array_like (float, 1D)", "The height matrix will be created in the same size as ``input``\n.. todo:: Explain, what height means in this case")
;
André Anjos's avatar
André Anjos committed
243

244
245
246
247
248
249
250
251
252
253
254
255
static xbob::extension::FunctionDoc s_norminv = xbob::extension::FunctionDoc(
  "norminv",
  "Computes the inverse normal cumulative distribution",
  "Computes the inverse normal cumulative distribution for a probability :math:`p`, given a distribution with mean :math:`\\mu` and standard deviation :math:`\\sigma`. "
  "Reference: http://home.online.no/~pjacklam/notes/invnorm/"
  )
  .add_prototype("p, mu, sigma", "inv")
  .add_parameter("p", "float", "The value to get the inverse distribution of, must lie in the range :math:`[0,1]`")
  .add_parameter("mu", "float", "The mean :math:`\\mu` of the normal distribution")
  .add_parameter("sigma", "float", "The standard deviation :math:`\\sigma` of the normal distribution")
  .add_return("inv", "float", "The inverse of the normal distribution")
;
André Anjos's avatar
André Anjos committed
256

257
258
259
260
261
262
263
264
265
266
267
static xbob::extension::FunctionDoc s_normsinv = xbob::extension::FunctionDoc(
  "normsinv",
  "Computes the inverse normal cumulative distribution",
  "Computes the inverse normal cumulative distribution for a probability :math:`p`, given a distribution with mean :math:`\\mu=0` and standard deviation :math:`\\sigma=1`. "
  "It is equivalent as calling ``norminv(p, 0, 1)`` (see :func:`norminv`). "
  "Reference: http://home.online.no/~pjacklam/notes/invnorm/"
  )
  .add_prototype("p", "inv")
  .add_parameter("p", "float", "The value to get the inverse distribution of, must lie in the range :math:`[0,1]`")
  .add_return("inv", "float", "The inverse of the normal distribution")
;
André Anjos's avatar
André Anjos committed
268

269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
static xbob::extension::FunctionDoc s_scatter = xbob::extension::FunctionDoc(
  "scatter",
  "Computes scatter matrix of a 2D array.",
  "Computes the scatter matrix of a 2D array *considering data is organized row-wise* (each sample is a row, each feature is a column). "
  "The resulting array ``s`` is squared with extents equal to the number of columns in ``a``. "
  "The resulting array ``m`` is a 1D array with the row means of ``a``. "
  "This function supports many calling modes, but you should provide, at least, the input data matrix ``a``. "
  "All non-provided arguments will be allocated internally and returned."
  )
  .add_prototype("a", "s, m")
  .add_prototype("a, s", "m")
  .add_prototype("a, m", "s")
  .add_prototype("a, s, m")
  .add_parameter("a", "array_like (float, 2D)", "The sample matrix, *considering data is organized row-wise* (each sample is a row, each feature is a column)")
  .add_parameter("s", "array_like (float, 2D)", "The scatter matrix, squared with extents equal to the number of columns in ``a``")
  .add_parameter("m", "array_like (float,1D)", "The mean matrix, with with the row means of ``a``")
  .add_return("s", "array_like (float, 2D)", "The scatter matrix, squared with extents equal to the number of columns in ``a``")
  .add_return("m", "array_like (float, 1D)", "The mean matrix, with with the row means of ``a``")
;
André Anjos's avatar
André Anjos committed
288

289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
static xbob::extension::FunctionDoc s_scatter_nocheck = xbob::extension::FunctionDoc(
  "scatter_",
  "Computes scatter matrix of a 2D array.",
  ".. warning:: This variant does not perform any checks on the input matrices and is faster then :func:`scatter`. "
  "Use it when you are sure your input matrices sizes match.\n"
  "Computes the scatter matrix of a 2D array *considering data is organized row-wise* (each sample is a row, each feature is a column). "
  "The resulting array ``s`` is squared with extents equal to the number of columns in ``a``. "
  "The resulting array ``m`` is a 1D array with the row means of ``a``. "
  "This function supports many calling modes, but you should provide, at least, the input data matrix ``a``. "
  "All non-provided arguments will be allocated internally and returned."
  )
  .add_prototype("a, s, m")
  .add_parameter("a", "array_like (float, 2D)", "The sample matrix, *considering data is organized row-wise* (each sample is a row, each feature is a column)")
  .add_parameter("s", "array_like (float, 2D)", "The scatter matrix, squared with extents equal to the number of columns in ``a``")
  .add_parameter("m", "array_like (float,1D)", "The mean matrix, with with the row means of ``a``")
;
André Anjos's avatar
André Anjos committed
305

306

307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
static xbob::extension::FunctionDoc s_scatters = xbob::extension::FunctionDoc(
  "scatters",
  "Computes :math:`S_w` and :math:`S_b` scatter matrices of a set of 2D arrays.",
  "Computes the within-class :math:`S_w` and between-class :math:`S_b` scatter matrices of a set of 2D arrays considering data is organized row-wise (each sample is a row, each feature is a column), and each matrix contains data of one class. "
  "Computes the scatter matrix of a 2D array *considering data is organized row-wise* (each sample is a row, each feature is a column). "
  "The implemented strategy is:\n"
  "1. Evaluate the overall mean (``m``), class means (:math:`m_k`) and the  total class counts (:math:`N`).\n"
  "2. Evaluate ``sw`` and ``sb`` using normal loops.\n"
  "Note that in this implementation, ``sw`` and ``sb`` will be normalized by N-1 (number of samples) and K (number of classes). "
  "This procedure makes the eigen values scaled by (N-1)/K, effectively increasing their values. "
  "The main motivation for this normalization are numerical precision concerns with the increasing number of samples causing a rather large :math:`S_w` matrix. "
  "A normalization strategy mitigates this problem. "
  "The eigen vectors will see no effect on this normalization as they are normalized in the euclidean sense (:math:`||a|| = 1`) so that does not change those.\n"
  "This function supports many calling modes, but you should provide, at least, the input ``data``. "
  "All non-provided arguments will be allocated internally and returned."
  )
  .add_prototype("data", "sw, sb, m")
  .add_prototype("data, sw, sb", "m")
  .add_prototype("data, sw, sb, m")
  .add_parameter("data", "[array_like (float, 2D)]", "The list of sample matrices. "
      "In each sample matrix the data is organized row-wise (each sample is a row, each feature is a column). "
      "Each matrix stores the data of a particular class. "
      "**Every matrix in ``data`` must have exactly the same number of columns.**")
  .add_parameter("sw", "array_like (float, 2D)", "The within-class scatter matrix :math:`S_w`, squared with extents equal to the number of columns in ``data``")
  .add_parameter("sb", "array_like (float, 2D)", "The between-class scatter matrix :math:`S_b`, squared with extents equal to the number of columns in ``data``")
  .add_parameter("m", "array_like (float,1D)", "The mean matrix, representing the ensemble mean with no prior (i.e., biased towards classes with more samples)")
  .add_return("sw", "array_like (float, 2D)", "The within-class scatter matrix :math:`S_w`")
  .add_return("sb", "array_like (float, 2D)", "The between-class scatter matrix :math:`S_b`")
  .add_return("m", "array_like (float, 1D)", "The mean matrix, representing the ensemble mean with no prior (i.e., biased towards classes with more samples)")
;
337

338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
static xbob::extension::FunctionDoc s_scatters_nocheck = xbob::extension::FunctionDoc(
  "scatters_",
  "Computes :math:`S_w` and :math:`S_b` scatter matrices of a set of 2D arrays.",
  ".. warning:: This variant does not perform any checks on the input matrices and is faster then :func:`scatters`. "
  "Use it when you are sure your input matrices sizes match.\n"
  "For a detailed description of the function, please see :func:`scatters`."
  )
  .add_prototype("data, sw, sb, m")
  .add_prototype("data, sw, sb")
  .add_parameter("data", "[array_like (float, 2D)]", "The list of sample matrices. "
      "In each sample matrix the data is organized row-wise (each sample is a row, each feature is a column). "
      "Each matrix stores the data of a particular class. "
      "**Every matrix in ``data`` must have exactly the same number of columns.**")
  .add_parameter("sw", "array_like (float, 2D)", "The within-class scatter matrix :math:`S_w`, squared with extents equal to the number of columns in ``data``")
  .add_parameter("sb", "array_like (float, 2D)", "The between-class scatter matrix :math:`S_b`, squared with extents equal to the number of columns in ``data``")
  .add_parameter("m", "array_like (float,1D)", "The mean matrix, representing the ensemble mean with no prior (i.e., biased towards classes with more samples)")
;
355
356


357
358
static PyMethodDef module_methods[] = {
    {
359
      s_histogram_intersection.name(),
360
361
      (PyCFunction)py_histogram_intersection,
      METH_VARARGS|METH_KEYWORDS,
362
      s_histogram_intersection.doc()
363
364
    },
    {
365
      s_chi_square.name(),
366
367
      (PyCFunction)py_chi_square,
      METH_VARARGS|METH_KEYWORDS,
368
      s_chi_square.doc()
369
370
    },
    {
371
      s_kullback_leibler.name(),
372
373
      (PyCFunction)py_kullback_leibler,
      METH_VARARGS|METH_KEYWORDS,
374
      s_kullback_leibler.doc()
375
    },
André Anjos's avatar
André Anjos committed
376
    {
377
      s_linsolve.name(),
André Anjos's avatar
André Anjos committed
378
379
      (PyCFunction)py_linsolve,
      METH_VARARGS|METH_KEYWORDS,
380
      s_linsolve.doc()
André Anjos's avatar
André Anjos committed
381
382
    },
    {
383
      s_linsolve_nocheck.name(),
André Anjos's avatar
André Anjos committed
384
385
      (PyCFunction)py_linsolve_nocheck,
      METH_VARARGS|METH_KEYWORDS,
386
      s_linsolve_nocheck.doc()
André Anjos's avatar
André Anjos committed
387
388
    },
    {
389
      s_linsolve_sympos.name(),
André Anjos's avatar
André Anjos committed
390
391
      (PyCFunction)py_linsolve_sympos,
      METH_VARARGS|METH_KEYWORDS,
392
      s_linsolve_sympos.doc()
André Anjos's avatar
André Anjos committed
393
394
    },
    {
395
      s_linsolve_sympos_nocheck.name(),
André Anjos's avatar
André Anjos committed
396
397
      (PyCFunction)py_linsolve_sympos_nocheck,
      METH_VARARGS|METH_KEYWORDS,
398
      s_linsolve_sympos_nocheck.doc()
André Anjos's avatar
André Anjos committed
399
400
    },
    {
401
      s_linsolve_cg_sympos.name(),
André Anjos's avatar
André Anjos committed
402
403
      (PyCFunction)py_linsolve_cg_sympos,
      METH_VARARGS|METH_KEYWORDS,
404
      s_linsolve_cg_sympos.doc()
André Anjos's avatar
André Anjos committed
405
406
    },
    {
407
      s_linsolve_cg_sympos_nocheck.name(),
André Anjos's avatar
André Anjos committed
408
409
      (PyCFunction)py_linsolve_cg_sympos_nocheck,
      METH_VARARGS|METH_KEYWORDS,
410
      s_linsolve_cg_sympos_nocheck.doc()
André Anjos's avatar
André Anjos committed
411
    },
André Anjos's avatar
André Anjos committed
412
    {
413
      s_pavx.name(),
André Anjos's avatar
André Anjos committed
414
415
      (PyCFunction)py_pavx,
      METH_VARARGS|METH_KEYWORDS,
416
      s_pavx.doc()
André Anjos's avatar
André Anjos committed
417
418
    },
    {
419
      s_pavx_nocheck.name(),
André Anjos's avatar
André Anjos committed
420
421
      (PyCFunction)py_pavx_nocheck,
      METH_VARARGS|METH_KEYWORDS,
422
      s_pavx_nocheck.doc()
André Anjos's avatar
André Anjos committed
423
424
    },
    {
425
      s_pavx_width.name(),
André Anjos's avatar
André Anjos committed
426
427
      (PyCFunction)py_pavx_width,
      METH_VARARGS|METH_KEYWORDS,
428
      s_pavx_width.doc()
André Anjos's avatar
André Anjos committed
429
430
    },
    {
431
      s_pavx_width_height.name(),
André Anjos's avatar
André Anjos committed
432
433
      (PyCFunction)py_pavx_width_height,
      METH_VARARGS|METH_KEYWORDS,
434
      s_pavx_width_height.doc()
André Anjos's avatar
André Anjos committed
435
    },
André Anjos's avatar
André Anjos committed
436
    {
437
      s_norminv.name(),
André Anjos's avatar
André Anjos committed
438
439
      (PyCFunction)py_norminv,
      METH_VARARGS|METH_KEYWORDS,
440
      s_norminv.doc()
André Anjos's avatar
André Anjos committed
441
442
    },
    {
443
      s_normsinv.name(),
André Anjos's avatar
André Anjos committed
444
445
      (PyCFunction)py_normsinv,
      METH_VARARGS|METH_KEYWORDS,
446
      s_normsinv.doc()
André Anjos's avatar
André Anjos committed
447
    },
448
    {
449
      s_scatter.name(),
450
451
      (PyCFunction)py_scatter,
      METH_VARARGS|METH_KEYWORDS,
452
      s_scatter.doc()
453
454
    },
    {
455
      s_scatter_nocheck.name(),
456
457
      (PyCFunction)py_scatter_nocheck,
      METH_VARARGS|METH_KEYWORDS,
458
      s_scatter_nocheck.doc()
459
460
    },
    {
461
      s_scatters.name(),
462
463
      (PyCFunction)py_scatters,
      METH_VARARGS|METH_KEYWORDS,
464
      s_scatters.doc()
465
466
    },
    {
467
      s_scatters_nocheck.name(),
468
469
      (PyCFunction)py_scatters_nocheck,
      METH_VARARGS|METH_KEYWORDS,
470
      s_scatters_nocheck.doc()
471
    },
472
473
474
475
476
    {0}  /* Sentinel */
};

PyDoc_STRVAR(module_docstr, "bob::math classes and methods");

André Anjos's avatar
André Anjos committed
477
478
479
480
481
482
#if PY_VERSION_HEX >= 0x03000000
static PyModuleDef module_definition = {
  PyModuleDef_HEAD_INIT,
  XBOB_EXT_MODULE_NAME,
  module_docstr,
  -1,
André Anjos's avatar
André Anjos committed
483
  module_methods,
André Anjos's avatar
André Anjos committed
484
485
486
487
  0, 0, 0, 0
};
#endif

André Anjos's avatar
André Anjos committed
488
static PyObject* create_module (void) {
489

490
  PyBobMathLpInteriorPoint_Type.tp_new = PyType_GenericNew;
André Anjos's avatar
André Anjos committed
491
  if (PyType_Ready(&PyBobMathLpInteriorPoint_Type) < 0) return 0;
492

493
  PyBobMathLpInteriorPointShortstep_Type.tp_base = &PyBobMathLpInteriorPoint_Type;
André Anjos's avatar
André Anjos committed
494
  if (PyType_Ready(&PyBobMathLpInteriorPointShortstep_Type) < 0) return 0;
495

496
  PyBobMathLpInteriorPointPredictorCorrector_Type.tp_base = &PyBobMathLpInteriorPoint_Type;
André Anjos's avatar
André Anjos committed
497
  if (PyType_Ready(&PyBobMathLpInteriorPointPredictorCorrector_Type) < 0) return 0;
498
499

  PyBobMathLpInteriorPointLongstep_Type.tp_base = &PyBobMathLpInteriorPoint_Type;
André Anjos's avatar
André Anjos committed
500
  if (PyType_Ready(&PyBobMathLpInteriorPointLongstep_Type) < 0) return 0;
André Anjos's avatar
André Anjos committed
501
502
503
504

# if PY_VERSION_HEX >= 0x03000000
  PyObject* m = PyModule_Create(&module_definition);
# else
André Anjos's avatar
André Anjos committed
505
  PyObject* m = Py_InitModule3(XBOB_EXT_MODULE_NAME, module_methods, module_docstr);
André Anjos's avatar
André Anjos committed
506
# endif
André Anjos's avatar
André Anjos committed
507
508
  if (!m) return 0;
  auto m_ = make_safe(m); ///< protects against early returns
509

André Anjos's avatar
André Anjos committed
510
511
512
  /* register version numbers and constants */
  if (PyModule_AddStringConstant(m, "__version__", XBOB_EXT_MODULE_VERSION) < 0)
    return 0;
513

514
515
  /* register the types to python */
  Py_INCREF(&PyBobMathLpInteriorPoint_Type);
André Anjos's avatar
André Anjos committed
516
517
  if (PyModule_AddObject(m, "LPInteriorPoint",
        (PyObject *)&PyBobMathLpInteriorPoint_Type) < 0) return 0;
518

519
  Py_INCREF(&PyBobMathLpInteriorPointShortstep_Type);
André Anjos's avatar
André Anjos committed
520
521
  if (PyModule_AddObject(m, "LPInteriorPointShortstep",
        (PyObject *)&PyBobMathLpInteriorPointShortstep_Type) < 0) return 0;
522

523
  Py_INCREF(&PyBobMathLpInteriorPointPredictorCorrector_Type);
André Anjos's avatar
André Anjos committed
524
525
  if (PyModule_AddObject(m, "LPInteriorPointPredictorCorrector",
        (PyObject *)&PyBobMathLpInteriorPointPredictorCorrector_Type) < 0) return 0;
526
527

  Py_INCREF(&PyBobMathLpInteriorPointLongstep_Type);
André Anjos's avatar
André Anjos committed
528
529
  if (PyModule_AddObject(m, "LPInteriorPointLongstep",
        (PyObject *)&PyBobMathLpInteriorPointLongstep_Type) < 0) return 0;
530
531

  /* imports xbob.blitz C-API */
André Anjos's avatar
André Anjos committed
532
  if (import_xbob_blitz() < 0) return 0;
533

André Anjos's avatar
André Anjos committed
534
  Py_INCREF(m);
André Anjos's avatar
André Anjos committed
535
536
  return m;

537
}
André Anjos's avatar
André Anjos committed
538
539
540
541
542
543
544

PyMODINIT_FUNC XBOB_EXT_ENTRY_NAME (void) {
# if PY_VERSION_HEX >= 0x03000000
  return
# endif
    create_module();
}