1 /*=================================================================================
2 * msvmocas_mex.c: OCAS solver for training multi-class linear SVM classifiers.
3 *
4 * Synopsis:
5 * [W,stat] = msvmocas(X,y,C,Method,TolRel,TolAbs,QPBound,BufSize,nExamples,MaxTime,verb)
6 *
7 * Input:
8 * X [nDim x nExamples] training feature inputs (sparse or dense matrix of doubles).
9 * y [nExamples x 1] labels; intgers 1,2,...nY
10 * C [1x1] regularization constant
11 * Method [1x1] 0..cutting plane; 1..OCAS (default 1)
12 * TolRel [1x1] halts if Q_P-Q_D <= abs(Q_P)*TolRel (default 0.01)
13 * TolAbs [1x1] halts if Q_P-Q_D <= TolAbs (default 0)
14 * QPValue [1x1] halts if Q_P <= QPBpound (default 0)
15 * BufSize [1x1] Initial size of active constrains buffer (default 2000)
16 * nExamples [1x1] Number of training examplesused for training; must be >0 and <= size(X,2).
17 * If nExamples = inf then nExamples is set to size(X,2).
18 * MaxTime [1x1] halts if time used by solver (data loading time is not counted) exceeds
19 * MaxTime given in seconds. Use MaxTime=inf (default) to switch off this stopping condition.
20 * verb [1x1] if non-zero then prints some info; (default 1)
21 *
22 * Output:
23 * W [nDim x nY] Paramater vectors of decision rule; [dummy,ypred] = max(W'*x)
24 * stat [struct] Optimizer statistics (field names are self-explaining).
25 *
26 * Copyright (C) 2008, 2012 Vojtech Franc, xfrancv@cmp.felk.cvut.cz
27 *
28 * This program is free software; you can redistribute it and/or
29 * modify it under the terms of the GNU General Public
30 * License as published by the Free Software Foundation;
31 *======================================================================================*/
32
33 #include <stdio.h>
34 #include <string.h>
35 #include <stdint.h>
36 #include <mex.h>
37
38 #include "libocas.h"
39 #include "ocas_helper.h"
40 #include "features_double.h"
41
42 #define DEFAULT_METHOD 1
43 #define DEFAULT_TOLREL 0.01
44 #define DEFAULT_TOLABS 0.0
45 #define DEFAULT_QPVALUE 0.0
46 #define DEFAULT_BUFSIZE 2000
47 #define DEFAULT_MAXTIME mxGetInf()
48 #define DEFAULT_VERB 1
49
50
51
52 /*======================================================================
53 Main code plus interface to Matlab.
54 ========================================================================*/
55
mexFunction(int nlhs,mxArray * plhs[],int nrhs,const mxArray * prhs[])56 void mexFunction( int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[] )
57 {
58 double C, TolRel, TolAbs, MaxTime, trn_err, QPBound;
59 double *ptr;
60 uint32_t i, j, BufSize;
61 uint16_t Method;
62 int verb;
63 ocas_return_value_T ocas;
64
65 /* timing variables */
66 double init_time;
67 double total_time;
68
69 total_time = get_time();
70 init_time = total_time;
71
72 if(nrhs < 3 || nrhs > 11)
73 mexErrMsgTxt("Improper number of input arguments.\n"
74 "\n"
75 "OCAS solver for training multi-class linear SVM classifiers.\n"
76 "\n"
77 "Synopsis:\n"
78 " [W,stat] = msvmocas(X,y,C,Method,TolRel,TolAbs,QPBound,BufSize,nExamples,MaxTime,verb)\n"
79 "\n"
80 "Input:\n"
81 " X [nDim x nExamples] training inputs (sparse or dense matrix of doubles).\n"
82 " y [nExamples x 1] labels must be integers 1,2,...nY\n"
83 " C [1x1] regularization constant\n"
84 " Method [1x1] 0..cutting plane; 1..OCAS (default 1)\n"
85 " TolRel [1x1] halts if Q_P-Q_D <= abs(Q_P)*TolRel (default 0.01)\n"
86 " TolAbs [1x1] halts if Q_P-Q_D <= TolAbs (default 0)\n"
87 " QPValue [1x1] halts if Q_P <= QPBpound (default 0)\n"
88 " BufSize [1x1] Initial size of active constrains buffer (default 2000)\n"
89 " nExamples [1x1] Number of training examples used for training; must be >0 and <= size(X,2).\n"
90 " If nExamples = inf then nExamples is set to size(X,2).\n"
91 " MaxTime [1x1] halts if time used by solver (data loading time is not counted) exceeds\n"
92 " MaxTime given in seconds. Use MaxTime=inf (default) to switch off this stopping condition.\n"
93 " verb [1x1] if non-zero then prints some info; (default 1)\n"
94 "\n"
95 "Output:\n"
96 " W [nDim x nY] Paramater vectors of decision rule; [dummy,ypred] = max(W'*x)\n"
97 " stat [struct] Optimizer statistics (field names are self-explaining).\n");
98
99 data_X = (mxArray*)prhs[0];
100 if( (mxGetNumberOfDimensions(data_X) != 2) ||
101 !( ( mxIsDouble(data_X) && mxIsSparse(data_X) ) ||
102 ( mxIsDouble(data_X) && !mxIsSparse(data_X) ) ))
103 {
104 mexErrMsgTxt("The first input argument must be two dimensional matrix of the following type:\n"
105 "dense double matrix or sparse double matrix.\n");
106 }
107
108 data_y = (double*)mxGetPr(prhs[1]);
109
110 if(LIBOCAS_MAX(mxGetM(prhs[1]),mxGetN(prhs[1])) != mxGetN(prhs[0]))
111 mexErrMsgTxt("Length of vector y must equal to the number of columns of matrix X.");
112
113 C = (double)mxGetScalar(prhs[2]);
114
115 if(nrhs >= 4)
116 Method = (uint32_t)mxGetScalar(prhs[3]);
117 else
118 Method = DEFAULT_METHOD;
119
120 if(nrhs >= 5)
121 TolRel = (double)mxGetScalar(prhs[4]);
122 else
123 TolRel = DEFAULT_TOLREL;
124
125 if(nrhs >= 6)
126 TolAbs = (double)mxGetScalar(prhs[5]);
127 else
128 TolAbs = DEFAULT_TOLABS;
129
130 if(nrhs >= 7)
131 QPBound = (double)mxGetScalar(prhs[6]);
132 else
133 QPBound = DEFAULT_QPVALUE;
134
135 if(nrhs >= 8)
136 BufSize = (uint32_t)mxGetScalar(prhs[7]);
137 else
138 BufSize = DEFAULT_BUFSIZE;
139
140 if(nrhs >= 9 && mxIsInf(mxGetScalar(prhs[8])) == false)
141 nData = (uint32_t)mxGetScalar(prhs[8]);
142 else
143 nData = mxGetN(data_X);
144
145 if(nData < 1 || nData > mxGetN(prhs[0]))
146 mexErrMsgTxt("Improper value of argument nData.");
147
148 if(nrhs >= 10)
149 MaxTime = (double)mxGetScalar(prhs[9]);
150 else
151 MaxTime = DEFAULT_MAXTIME;
152
153 if(nrhs >= 11)
154 verb = (int)mxGetScalar(prhs[10]);
155 else
156 verb = DEFAULT_VERB;
157
158
159 nDim = mxGetM(data_X);
160 for(i=0, nY = 0; i < nData; i++)
161 nY = LIBOCAS_MAX(nY, (uint32_t)data_y[i]);
162
163 /*----------------------------------------------------------------
164 Print setting
165 -------------------------------------------------------------------*/
166 if(verb)
167 {
168 mexPrintf("Input data statistics:\n"
169 " # of examples : %d\n"
170 " # of classes : %d\n"
171 " dimensionality : %d\n",
172 nData, nY, nDim);
173
174 if( mxIsSparse(data_X)== true )
175 mexPrintf(" density : %.2f%%\n",
176 100.0*(double)mxGetNzmax(data_X)/((double)nDim*(double)(mxGetN(data_X))));
177 else
178 mexPrintf(" density : 100%% (full)\n");
179
180 mexPrintf("Setting:\n"
181 " C : %f\n"
182 " # of examples : %d\n"
183 " solver : %d\n"
184 " cache size : %d\n"
185 " TolAbs : %f\n"
186 " TolRel : %f\n"
187 " QPValue : %f\n"
188 " MaxTime : %f [s]\n",
189 C, nData, Method,BufSize,TolAbs,TolRel, QPBound, MaxTime);
190 }
191
192 /* learned weight vector */
193 plhs[0] = (mxArray*)mxCreateDoubleMatrix(nDim,nY,mxREAL);
194 W = (double*)mxGetPr(plhs[0]);
195 if(W == NULL) mexErrMsgTxt("Not enough memory for vector W.");
196
197 oldW = (double*)mxCalloc(nY*nDim,sizeof(double));
198 if(oldW == NULL) mexErrMsgTxt("Not enough memory for vector oldW.");
199
200 /* allocate buffer for computing cutting plane */
201 new_a = (double*)mxCalloc(nY*nDim,sizeof(double));
202 if(new_a == NULL)
203 mexErrMsgTxt("Not enough memory for auxciliary cutting plane buffer new_a.");
204
205 /* select function to print progress info */
206 void (*print_function)(ocas_return_value_T);
207 if(verb)
208 {
209 mexPrintf("Starting optimization:\n");
210 print_function = &ocas_print;
211 }
212 else
213 {
214 print_function = &ocas_print_null;
215 }
216
217
218 if( mxIsSparse(data_X)== true )
219 {
220 /* init cutting plane buffer */
221 sparse_A.nz_dims = mxCalloc(BufSize,sizeof(uint32_t));
222 sparse_A.index = mxCalloc(BufSize,sizeof(sparse_A.index[0]));
223 sparse_A.value = mxCalloc(BufSize,sizeof(sparse_A.value[0]));
224 if(sparse_A.nz_dims == NULL || sparse_A.index == NULL || sparse_A.value == NULL)
225 mexErrMsgTxt("Not enough memory for cutting plane buffer sparse_A.");
226
227 init_time=get_time()-init_time;
228
229 ocas = msvm_ocas_solver( C, data_y, nY, nData, TolRel, TolAbs, QPBound, MaxTime,BufSize, Method,
230 &msvm_sparse_compute_W, &msvm_update_W, &msvm_sparse_add_new_cut,
231 &msvm_sparse_compute_output, &qsort_data, print_function, 0);
232 }
233 else
234 {
235 /* init cutting plane buffer */
236 full_A = mxCalloc(BufSize*nDim*nY,sizeof(double));
237 if( full_A == NULL )
238 mexErrMsgTxt("Not enough memory for cutting plane buffer full_A.");
239
240 init_time=get_time()-init_time;
241
242 ocas = msvm_ocas_solver( C, data_y, nY, nData, TolRel, TolAbs, QPBound, MaxTime,BufSize, Method,
243 &msvm_full_compute_W, &msvm_update_W, &msvm_full_add_new_cut,
244 &msvm_full_compute_output, &qsort_data, print_function, 0);
245 }
246
247 total_time=get_time()-total_time;
248
249 if(verb)
250 {
251 mexPrintf("Stopping condition: ");
252 switch( ocas.exitflag )
253 {
254 case 1: mexPrintf("1-Q_D/Q_P <= TolRel(=%f) satisfied.\n", TolRel); break;
255 case 2: mexPrintf("Q_P-Q_D <= TolAbs(=%f) satisfied.\n", TolAbs); break;
256 case 3: mexPrintf("Q_P <= QPBound(=%f) satisfied.\n", QPBound); break;
257 case 4: mexPrintf("Optimization time (=%f) >= MaxTime(=%f).\n", ocas.ocas_time, MaxTime); break;
258 case -1: mexPrintf("Has not converged!\n" ); break;
259 case -2: mexPrintf("Not enough memory for the solver.\n" ); break;
260 }
261
262 mexPrintf("Timing statistics:\n"
263 " init_time : %f[s]\n"
264 " qp_solver_time : %f[s]\n"
265 " sort_time : %f[s]\n"
266 " output_time : %f[s]\n"
267 " add_time : %f[s]\n"
268 " w_time : %f[s]\n"
269 " print_time : %f[s]\n"
270 " ocas_time : %f[s]\n"
271 " total_time : %f[s]\n",
272 init_time, ocas.qp_solver_time, ocas.sort_time, ocas.output_time,
273 ocas.add_time, ocas.w_time, ocas.print_time, ocas.ocas_time, total_time);
274
275 mexPrintf("Training error: %.4f%%\n", 100*(double)ocas.trn_err/(double)nData);
276 }
277
278 const char *field_names[] = {"nTrnErrors","Q_P","Q_D","nIter","nCutPlanes","exitflag",
279 "init_time","output_time","sort_time","qp_solver_time","add_time",
280 "w_time","ocas_time","total_time"};
281 mwSize dims[2] = {1,1};
282
283 plhs[1] = mxCreateStructArray(2, dims, (sizeof(field_names)/sizeof(*field_names)), field_names);
284
285 mxSetField(plhs[1],0,"nIter",mxCreateDoubleScalar((double)ocas.nIter));
286 mxSetField(plhs[1],0,"nCutPlanes",mxCreateDoubleScalar((double)ocas.nCutPlanes));
287 mxSetField(plhs[1],0,"nTrnErrors",mxCreateDoubleScalar(ocas.trn_err));
288 mxSetField(plhs[1],0,"Q_P",mxCreateDoubleScalar(ocas.Q_P));
289 mxSetField(plhs[1],0,"Q_D",mxCreateDoubleScalar(ocas.Q_D));
290 mxSetField(plhs[1],0,"init_time",mxCreateDoubleScalar(init_time));
291 mxSetField(plhs[1],0,"output_time",mxCreateDoubleScalar(ocas.output_time));
292 mxSetField(plhs[1],0,"sort_time",mxCreateDoubleScalar(ocas.sort_time));
293 mxSetField(plhs[1],0,"qp_solver_time",mxCreateDoubleScalar(ocas.qp_solver_time));
294 mxSetField(plhs[1],0,"add_time",mxCreateDoubleScalar(ocas.add_time));
295 mxSetField(plhs[1],0,"w_time",mxCreateDoubleScalar(ocas.w_time));
296 mxSetField(plhs[1],0,"ocas_time",mxCreateDoubleScalar(ocas.ocas_time));
297 mxSetField(plhs[1],0,"total_time",mxCreateDoubleScalar(total_time));
298 mxSetField(plhs[1],0,"exitflag",mxCreateDoubleScalar((double)ocas.exitflag));
299
300 return;
301 }
302
303