[libocas] 55/60: Imported Upstream version 0.97
Christian Kastner
chrisk-guest at moszumanska.debian.org
Mon Aug 25 03:34:46 UTC 2014
This is an automated email from the git hooks/post-receive script.
chrisk-guest pushed a commit to branch master
in repository libocas.
commit 06eb50d88fa2b22cd4b0055e16b8b23583cbd746
Author: Christian Kastner <debian at kvr.at>
Date: Wed Feb 19 15:43:36 2014 +0100
Imported Upstream version 0.97
---
ChangeLog | 34 +-
Contents.m | 1 +
Makefile | 26 +-
RELEASES | 3 -
compute_auc.mexa64 | Bin 27573 -> 0 bytes
features_bool.c | 177 +++++++++++
features_bool.h | 20 ++
features_single.c | 2 +-
html/ChangeLog | 27 +-
html/index.html | 11 +-
lbppyr_features.mexa64 | Bin 25711 -> 0 bytes
libocas.c | 492 ++++++++++++++++++++++++++++-
libocas.h | 24 ++
libocas.so | Bin 41123 -> 0 bytes
libqp_splx.c | 814 ++++++++++++++++++++++++------------------------
linclassif | Bin 17743 -> 0 bytes
linclassif_light.mexa64 | Bin 31853 -> 0 bytes
load_svmlight_file.c | 244 +++++++++++++++
msvmocas | Bin 78218 -> 0 bytes
msvmocas.mexa64 | Bin 65734 -> 0 bytes
msvmocas_light.mexa64 | Bin 65780 -> 0 bytes
ocas_helper.c | 47 +++
ocas_helper.h | 10 +-
svmocas | Bin 82262 -> 0 bytes
svmocas.mexa64 | Bin 70189 -> 0 bytes
svmocas_bool_mex.c | 314 +++++++++++++++++++
svmocas_bool_test.m | 16 +
svmocas_lbp.mexa64 | Bin 51481 -> 0 bytes
svmocas_light.mexa64 | Bin 65826 -> 0 bytes
svmocas_nnw.m | 74 +++++
svmocas_nnw_mex.c | 456 +++++++++++++++++++++++++++
version.h | 2 +-
32 files changed, 2339 insertions(+), 455 deletions(-)
diff --git a/ChangeLog b/ChangeLog
index f20ff02..dcf8e37 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,17 +1,27 @@
+2013-05-14 Vojtech Franc
+ * Bug in Makefile fixed (reported by fhy881229 at zju.edu.cn)
+ * Release v0.97 (SVN revision 42)
+ * Added Matlab function for loading data in SVM light format.
+2012-07-23 Vojtech Franc
+ * Release of v0.95
+ * Added solver SVMOCAS_NNW which allows training two-cclass SVM
+ * with additinal constraints enforsing a selected subset of weights
+ * to be non-negative. There is currently only a Matlab MEX interface for
+ * this solver.
2012-02-06 Vojtech Franc
- * Release of v0.95
- * Helper functions for different data types put to separate files -
- * features_double.c/h, features_single.c/h, features_int8.c.h.
- * Programmed support for data stored in int8 and float single precision.
- * Separate MEX funtions for training from SVM^light format.
- * Removed unused functions for computing LBP features.
+ * Release of v0.95
+ * Helper functions for different data types put to separate files -
+ * features_double.c/h, features_single.c/h, features_int8.c.h.
+ * Programmed support for data stored in int8 and float single precision.
+ * Separate MEX funtions for training from SVM^light format.
+ * Removed unused functions for computing LBP features.
2010-06-11 Vojtech Franc
- * Added functions which implement the COFFIN framework for training
- * translation invariant image classifier from virtual example.
- * Added functions:
- * lbppyr.m, lbppyr_mex.c, liblbp.h, ocas_lbp_helper.h, svmocas_lbp.m
- * lbpfilter_mex.c, lbppyr_features_mex.c, liblbp.c, ocas_lbp_helper.c
- * svmocas_lbp_example.m, svmocas_lbp_mex.c
+ * Added functions which implement the COFFIN framework for training
+ * translation invariant image classifier from virtual example.
+ * Added functions:
+ * lbppyr.m, lbppyr_mex.c, liblbp.h, ocas_lbp_helper.h, svmocas_lbp.m
+ * lbpfilter_mex.c, lbppyr_features_mex.c, liblbp.c, ocas_lbp_helper.c
+ * svmocas_lbp_example.m, svmocas_lbp_mex.c
2009-08-03 Vojtech Franc
* BUG FIX: OCAS solver was crashing on some data (reported by Alex Binder).
The problem was a bug in "qsort_data" function which did not expect a sequnce of length 1
diff --git a/Contents.m b/Contents.m
index deb8ec7..d1be82a 100644
--- a/Contents.m
+++ b/Contents.m
@@ -7,6 +7,7 @@
% SVM solvers for training linear two-class classifiers:
% svmocas Accepts examples stored in dense double or sparse
% double or dense single or dense int8 matrix.
+% svmocas_nnw Allows additional constrains enforcing non-negative weights.
% svmocas_light Loads examples from file in SVM^light format.
% svmocas_lbp Examples are LBP features computed on a set of
% given grayscale images.
diff --git a/Makefile b/Makefile
index 6bc4c9a..f847bf9 100644
--- a/Makefile
+++ b/Makefile
@@ -11,13 +11,16 @@ CLIBS := -lm
ifeq (yes,$(MEXDETECTED))
-all: svmocas.$(MEXSUFFIX) svmocas_light.$(MEXSUFFIX) linclassif_light.$(MEXSUFFIX) libocas.so svmocas msvmocas linclassif msvmocas.$(MEXSUFFIX) msvmocas_light.$(MEXSUFFIX) compute_auc.$(MEXSUFFIX) svmocas_lbp.$(MEXSUFFIX) lbppyr_features.$(MEXSUFFIX)
+all: svmocas_nnw.$(MEXSUFFIX) svmocas.$(MEXSUFFIX) svmocas_light.$(MEXSUFFIX) linclassif_light.$(MEXSUFFIX) libocas.so svmocas msvmocas linclassif msvmocas.$(MEXSUFFIX) msvmocas_light.$(MEXSUFFIX) compute_auc.$(MEXSUFFIX) svmocas_lbp.$(MEXSUFFIX) svmocas_bool.$(MEXSUFFIX) lbppyr_features.$(MEXSUFFIX) load_svmlight_file.$(MEXSUFFIX)
-compute_auc.$(MEXSUFFIX): compute_auc_mex.c libocas.h ocas_helper.h ocas_helper.c
+compute_auc.$(MEXSUFFIX): compute_auc_mex.c libocas.h ocas_helper.h ocas_helper.c version.h
$(MEX) $(MEXFLAGS) -DLIBOCAS_MATLAB -O -output compute_auc.$(MEXSUFFIX) compute_auc_mex.c ocas_helper.c lib_svmlight_format.c
svmocas.$(MEXSUFFIX): libocas.c libocas.h libqp_splx.c libqp.h svmocas_mex.c lib_svmlight_format.c lib_svmlight_format.h ocas_helper.c ocas_helper.h features_int8.h features_int8.c features_double.h features_double.c features_single.h features_single.c
- $(MEX) $(MEXFLAGS) -DLIBOCAS_MATLAB -O -output svmocas.$(MEXSUFFIX) svmocas_mex.c lib_svmlight_format.c ocas_helper.c features_int8.c features_double.c features_single.c libocas.c libqp_splx.c
+ $(MEX) $(MEXFLAGS) -DLIBOCAS_MATLAB -O -output svmocas.$(MEXSUFFIX) svmocas_mex.c ocas_helper.c lib_svmlight_format.c features_int8.c features_double.c features_single.c libocas.c libqp_splx.c
+
+svmocas_nnw.$(MEXSUFFIX): libocas.c libocas.h libqp_splx.c libqp.h svmocas_nnw_mex.c lib_svmlight_format.c lib_svmlight_format.h ocas_helper.c ocas_helper.h features_int8.h features_int8.c features_double.h features_double.c features_single.h features_single.c
+ $(MEX) $(MEXFLAGS) -DLIBOCAS_MATLAB -O -output svmocas_nnw.$(MEXSUFFIX) svmocas_nnw_mex.c lib_svmlight_format.c ocas_helper.c features_int8.c features_double.c features_single.c libocas.c libqp_splx.c
svmocas_light.$(MEXSUFFIX): libocas.c libocas.h libqp_splx.c libqp.h svmocas_light_mex.c lib_svmlight_format.c lib_svmlight_format.h ocas_helper.c ocas_helper.h features_double.h features_double.c
$(MEX) $(MEXFLAGS) -DLIBOCAS_MATLAB -O -output svmocas_light.$(MEXSUFFIX) svmocas_light_mex.c lib_svmlight_format.c ocas_helper.c features_double.c libocas.c libqp_splx.c
@@ -37,29 +40,34 @@ lbppyr.$(MEXSUFFIX): lbppyr_mex.c
svmocas_lbp.$(MEXSUFFIX): libocas.c libocas.h libqp_splx.c libqp.h svmocas_lbp_mex.c ocas_lbp_helper.c ocas_lbp_helper.h liblbp.h liblbp.c
$(MEX) $(MEXFLAGS) -DLIBOCAS_MATLAB -O -output svmocas_lbp.$(MEXSUFFIX) svmocas_lbp_mex.c ocas_lbp_helper.c libocas.c libqp_splx.c liblbp.c
+svmocas_bool.$(MEXSUFFIX): libocas.c libocas.h libqp_splx.c libqp.h svmocas_bool_mex.c features_bool.c features_bool.h liblbp.h liblbp.c
+ $(MEX) $(MEXFLAGS) -DLIBOCAS_MATLAB -O -output svmocas_bool.$(MEXSUFFIX) svmocas_bool_mex.c features_bool.c ocas_helper.c lib_svmlight_format.c libocas.c libqp_splx.c
+
lbppyr_features.$(MEXSUFFIX): lbppyr_features_mex.c liblbp.c liblbp.h
$(MEX) -g $(MEXFLAGS) -DLIBOCAS_MATLAB -O -output lbppyr_features.$(MEXSUFFIX) lbppyr_features_mex.c liblbp.c
lbpfilter.$(MEXSUFFIX): lbpfilter_mex.c
$(MEX) $(MEXFLAGS) -DLIBOCAS_MATLAB -O -output lbpfilter.$(MEXSUFFIX) lbpfilter_mex.c
+load_svmlight_file.$(MEXSUFFIX): load_svmlight_file.c lib_svmlight_format.c lib_svmlight_format.h
+ $(MEX) $(MEXFLAGS) -DLIBOCAS_MATLAB -O -output load_svmlight_file.$(MEXSUFFIX) load_svmlight_file.c lib_svmlight_format.c
else
-all: libocas.so svmocas msvmocas linclass
+all: libocas.so svmocas msvmocas linclassif
endif
-svmocas: svmocas.c lib_svmlight_format.c sparse_mat.c ocas_helper.c ocas_helper.h libocas.h sparse_mat.h libocas.c features_double.h features_double.c
+svmocas: svmocas.c svmocas.c lib_svmlight_format.c sparse_mat.c ocas_helper.c ocas_helper.h libocas.h sparse_mat.h libocas.c features_double.h features_double.c version.h
$(CC) $(CFLAGS) -o $@ svmocas.c lib_svmlight_format.c sparse_mat.c ocas_helper.c features_double.c libocas.c libqp_splx.c $(CLIBS)
-msvmocas: msvmocas.c lib_svmlight_format.c sparse_mat.c ocas_helper.c ocas_helper.h libocas.h sparse_mat.h libocas.c features_double.h features_double.c
+msvmocas: msvmocas.c lib_svmlight_format.c sparse_mat.c ocas_helper.c ocas_helper.h libocas.h sparse_mat.h libocas.c features_double.h features_double.c version.h
$(CC) $(CFLAGS) -o $@ msvmocas.c lib_svmlight_format.c sparse_mat.c ocas_helper.c features_double.c libocas.c libqp_splx.c $(CLIBS)
-linclassif: linclassif.c lib_svmlight_format.c libocas.h
+linclassif: linclassif.c lib_svmlight_format.c libocas.h version.h
$(CC) $(CFLAGS) -o $@ linclassif.c lib_svmlight_format.c $(CLIBS)
-libocas.so: libocas.c libocas.h libqp_splx.c libqp.h
+libocas.so: libocas.c libocas.h libqp_splx.c libqp.h
$(CC) $(CFLAGS) -shared -o $@ libocas.c libqp_splx.c $(CLIBS)
clean:
- rm -f *~ svmocas.$(MEXSUFFIX) svmocas_light.$(MEXSUFFIX) linclassif_light.$(MEXSUFFIX) svmocas msvmocas linclassif libocas.so msvmocas.$(MEXSUFFIX) msvmocas_light.$(MEXSUFFIX) lbppyr_features.$(MEXSUFFIX) svmocas_lbp.$(MEXSUFFIX) compute_auc.$(MEXSUFFIX)
+ rm -f *~ svmocas.$(MEXSUFFIX) svmocas_nnw.$(MEXSUFFIX) svmocas_bool.$(MEXSUFFIX) svmocas_light.$(MEXSUFFIX) linclassif_light.$(MEXSUFFIX) svmocas msvmocas linclassif libocas.so msvmocas.$(MEXSUFFIX) msvmocas_light.$(MEXSUFFIX) lbppyr_features.$(MEXSUFFIX) svmocas_lbp.$(MEXSUFFIX) compute_auc.$(MEXSUFFIX) load_svmlight_file.$(MEXSUFFIX)
diff --git a/RELEASES b/RELEASES
deleted file mode 100644
index d7508b0..0000000
--- a/RELEASES
+++ /dev/null
@@ -1,3 +0,0 @@
-Version 0.92 2009-08-03 SVN Revision 770
-Version 0.91 2009-07-14 SVN Revision 760
-Version 0.9 2008-12-10 SVN Revision 654
diff --git a/compute_auc.mexa64 b/compute_auc.mexa64
deleted file mode 100755
index 975441d..0000000
Binary files a/compute_auc.mexa64 and /dev/null differ
diff --git a/features_bool.c b/features_bool.c
new file mode 100644
index 0000000..ff7d754
--- /dev/null
+++ b/features_bool.c
@@ -0,0 +1,177 @@
+/*-----------------------------------------------------------------------
+ * ocas_helper.c: Implementation of helper functions for the OCAS solver.
+ *
+ * It supports both sparse and dense matrices and loading data from
+ * the SVM^light format.
+ *-------------------------------------------------------------------- */
+
+#define _FILE_OFFSET_BITS 64
+
+#include <pthread.h>
+
+#include <stdio.h>
+#include <string.h>
+#include <stdint.h>
+#include <sys/time.h>
+#include <stdlib.h>
+#include <time.h>
+
+#include "libocas.h"
+#include "ocas_helper.h"
+
+
+/*----------------------------------------------------------------------
+ full_bool_compute_output( output ) does the follwing:
+
+ output = data_X'*W;
+ ----------------------------------------------------------------------*/
+int full_bool_compute_output( double *output, void* user_data )
+{
+ uint32_t i, j, k, dataDim; /* k is new index for the vector W */
+ double tmp;
+ uint8_t* ptr;
+
+ ptr = (uint8_t*)mxGetPr( data_X );
+
+ dataDim = nDim/8;
+
+ for(i=0; i < nData; i++) {
+ tmp = X0*W0;
+
+ k = 0;
+
+ for(j=0; j < dataDim; j++ ) {
+ if (ptr[LIBOCAS_INDEX(j,i,dataDim)] & 0x01) {
+ tmp += W[k];
+ }
+ k++;
+ if (ptr[LIBOCAS_INDEX(j,i,dataDim)] & 0x02) {
+ tmp += W[k];
+ }
+ k++;
+ if (ptr[LIBOCAS_INDEX(j,i,dataDim)] & 0x04) {
+ tmp += W[k];
+ }
+ k++;
+ if (ptr[LIBOCAS_INDEX(j,i,dataDim)] & 0x08) {
+ tmp += W[k];
+ }
+ k++;
+ if (ptr[LIBOCAS_INDEX(j,i,dataDim)] & 0x10) {
+ tmp += W[k];
+ }
+ k++;
+ if (ptr[LIBOCAS_INDEX(j,i,dataDim)] & 0x20) {
+ tmp += W[k];
+ }
+ k++;
+ if (ptr[LIBOCAS_INDEX(j,i,dataDim)] & 0x40) {
+ tmp += W[k];
+ }
+ k++;
+ if (ptr[LIBOCAS_INDEX(j,i,dataDim)] & 0x80) {
+ tmp += W[k];
+ }
+ k++;
+ }
+ output[i] = tmp*data_y[i];
+ }
+
+ return 0;
+}
+
+
+
+/*----------------------------------------------------------------------------------
+ full_bool_add_new_cut( new_col_H, new_cut, cut_length, nSel ) does the following:
+
+ new_a = sum(data_X(:,find(new_cut ~=0 )),2);
+ new_col_H = [full_A(:,1:nSel)'*new_a ; new_a'*new_a];
+ full_A(:,nSel+1) = new_a;
+
+ ---------------------------------------------------------------------------------*/
+int full_bool_add_new_cut( double *new_col_H,
+ uint32_t *new_cut,
+ uint32_t cut_length,
+ uint32_t nSel,
+ void* user_data)
+{
+ /*****************************/
+ double sq_norm_a;
+ uint32_t i, j, k, dataDim; /* k is new index for the vector a*/
+ uint8_t *ptr;
+
+ ptr = (uint8_t*)mxGetPr(data_X);
+
+ dataDim = nDim/8;
+
+ memset(new_a, 0, sizeof(double)*nDim);
+
+
+ for(i=0; i < cut_length; i++) {
+
+ k = 0;
+
+ for(j=0; j < dataDim; j++ ) {
+ if (ptr[LIBOCAS_INDEX(j,new_cut[i],dataDim)] & 0x01) {
+ new_a[k] += data_y[new_cut[i]];
+ }
+ k++;
+ if (ptr[LIBOCAS_INDEX(j,new_cut[i],dataDim)] & 0x02) {
+ new_a[k] += data_y[new_cut[i]];
+ }
+ k++;
+ if (ptr[LIBOCAS_INDEX(j,new_cut[i],dataDim)] & 0x04) {
+ new_a[k] += data_y[new_cut[i]];
+ }
+ k++;
+ if (ptr[LIBOCAS_INDEX(j,new_cut[i],dataDim)] & 0x08) {
+ new_a[k] += data_y[new_cut[i]];
+ }
+ k++;
+ if (ptr[LIBOCAS_INDEX(j,new_cut[i],dataDim)] & 0x10) {
+ new_a[k] += data_y[new_cut[i]];
+ }
+ k++;
+ if (ptr[LIBOCAS_INDEX(j,new_cut[i],dataDim)] & 0x20) {
+ new_a[k] += data_y[new_cut[i]];
+ }
+ k++;
+ if (ptr[LIBOCAS_INDEX(j,new_cut[i],dataDim)] & 0x40) {
+ new_a[k] += data_y[new_cut[i]];
+ }
+ k++;
+ if (ptr[LIBOCAS_INDEX(j,new_cut[i],dataDim)] & 0x80) {
+ new_a[k] += data_y[new_cut[i]];
+ }
+ k++;
+ }
+
+ A0[nSel] += X0*data_y[new_cut[i]];
+ }
+
+ /*****************************/
+
+ /* compute new_a'*new_a and insert new_a to the last column of full_A */
+ sq_norm_a = A0[nSel]*A0[nSel];
+ for(j=0; j < nDim; j++ ) {
+ sq_norm_a += new_a[j]*new_a[j];
+ full_A[LIBOCAS_INDEX(j,nSel,nDim)] = new_a[j];
+ }
+
+ new_col_H[nSel] = sq_norm_a;
+ for(i=0; i < nSel; i++) {
+ double tmp = A0[nSel]*A0[i];
+
+ for(j=0; j < nDim; j++ ) {
+ tmp += new_a[j]*full_A[LIBOCAS_INDEX(j,i,nDim)];
+ }
+ new_col_H[i] = tmp;
+ }
+
+/* mxFree( new_a );*/
+
+ return 0;
+}
+
+
diff --git a/features_bool.h b/features_bool.h
new file mode 100644
index 0000000..579907e
--- /dev/null
+++ b/features_bool.h
@@ -0,0 +1,20 @@
+/*-----------------------------------------------------------------------
+ * ocas_helper.h: Implementation of helper functions for the OCAS solver.
+ *
+ * It supports both sparse and dense matrices and loading data from
+ * the SVM^light format.
+ *-------------------------------------------------------------------- */
+
+#ifndef _svmocas_bool_helper_h
+#define _svmocas_bool_helper_h
+
+
+int full_bool_compute_output( double *output, void* user_data );
+int full_bool_add_new_cut( double *new_col_H,
+ uint32_t *new_cut,
+ uint32_t cut_length,
+ uint32_t nSel,
+ void* user_data);
+
+
+#endif
diff --git a/features_single.c b/features_single.c
index 0b8296e..4117235 100644
--- a/features_single.c
+++ b/features_single.c
@@ -88,4 +88,4 @@ int full_single_compute_output( double *output, void* user_data )
}
return 0;
-}
\ No newline at end of file
+}
diff --git a/html/ChangeLog b/html/ChangeLog
index cd5b4d2..c66654d 100644
--- a/html/ChangeLog
+++ b/html/ChangeLog
@@ -1,12 +1,23 @@
-2011-12-22 Vojtech Fracn
- * Added "svmocas_int8" supporting features represented as int8.
+2012-07-23 Vojtech Franc
+ * Release of v0.95
+ * Added solver SVMOCAS_NNW which allows training two-cclass SVM
+ * with additinal constraints enforsing a selected subset of weights
+ * to be non-negative. There is currently only a Matlab MEX interface for
+ * this solver.
+2012-02-06 Vojtech Franc
+ * Release of v0.95
+ * Helper functions for different data types put to separate files -
+ * features_double.c/h, features_single.c/h, features_int8.c.h.
+ * Programmed support for data stored in int8 and float single precision.
+ * Separate MEX funtions for training from SVM^light format.
+ * Removed unused functions for computing LBP features.
2010-06-11 Vojtech Franc
- * Added functions which implement the COFFIN framework for training
- * translation invariant image classifier from virtual example.
- * Added functions:
- * lbppyr.m, lbppyr_mex.c, liblbp.h, ocas_lbp_helper.h, svmocas_lbp.m
- * lbpfilter_mex.c, lbppyr_features_mex.c, liblbp.c, ocas_lbp_helper.c
- * svmocas_lbp_example.m, svmocas_lbp_mex.c
+ * Added functions which implement the COFFIN framework for training
+ * translation invariant image classifier from virtual example.
+ * Added functions:
+ * lbppyr.m, lbppyr_mex.c, liblbp.h, ocas_lbp_helper.h, svmocas_lbp.m
+ * lbpfilter_mex.c, lbppyr_features_mex.c, liblbp.c, ocas_lbp_helper.c
+ * svmocas_lbp_example.m, svmocas_lbp_mex.c
2009-08-03 Vojtech Franc
* BUG FIX: OCAS solver was crashing on some data (reported by Alex Binder).
The problem was a bug in "qsort_data" function which did not expect a sequnce of length 1
diff --git a/html/index.html b/html/index.html
index 4d8814c..a292d38 100644
--- a/html/index.html
+++ b/html/index.html
@@ -13,7 +13,7 @@ classifiers from large-scale data</H2>
<A href="http://ida.first.fraunhofer.de/homepages/sonne/first/">Soeren
Sonnenburg</a>
<p>
-Last Modified: 06-Feb-2012
+Last Modified: 23-Jul-2012
</div>
@@ -51,6 +51,7 @@ notebook with just 4GB of memory.
<li>Two-class and genuine multi-class SVM formulations.</li>
<li>A stand alone application and MEX interface for Matlab.</li>
<li>Optimized for features represented by different data types: sparse/dense, double/single precision, int8.</li>
+ <li>Allows training two-class SVM classifier with additional constraints enforcing a subset of selected weights to be non-negative.</li>
<li>Reads examples from SVM^light format.</li>
<li>Parallelized version of the two-class SVM solver.</li>
<li>Allows using different C for each training example (only for Matlab).</li>
@@ -90,7 +91,13 @@ from <img align="top" border="0" src="set_y_sc.png">.
LIBOCAS can be downloaded from here:
<ul>
- <li><img align="top" border="0" src="new.gif"> Version 0.95, 2012-02-06, <a
+ <li><img align="top" border="0" src="new.gif"> Version 0.97, 2013-05-14, <a
+ href="http://cmp.felk.cvut.cz/~xfrancv/ocas/libocas_v097.zip">libocas_v097.zip</a>: Bug in makefile fixed. </li>
+ <!-- ------>
+ <li>Version 0.96, 2012-07-23, <a
+ href="http://cmp.felk.cvut.cz/~xfrancv/ocas/libocas_v096.zip">libocas_v096.zip</a>: New solver which enforses a subset of selected weights to be non-negative
+ </li>
+ <li>Version 0.95, 2012-02-06, <a
href="http://cmp.felk.cvut.cz/~xfrancv/ocas/libocas_v095.zip">libocas_v095.zip</a>: Support for features in single precision floats. Code polished. </li>
<li>Version 0.94, 2011-12-22, <a
href="http://cmp.felk.cvut.cz/~xfrancv/ocas/libocas_v094.zip">libocas_v094.zip</a>: Support for features represented as int8. </li>
diff --git a/lbppyr_features.mexa64 b/lbppyr_features.mexa64
deleted file mode 100755
index 65c4668..0000000
Binary files a/lbppyr_features.mexa64 and /dev/null differ
diff --git a/libocas.c b/libocas.c
index b1d0182..2e6757b 100644
--- a/libocas.c
+++ b/libocas.c
@@ -22,7 +22,8 @@
#include "libocas.h"
#include "libqp.h"
-#define LAMBDA 0.1 /* must be from (0,1> 1..means that OCAS becomes equivalent to CPA */
+#define MU 0.1 /* must be from (0,1> 1..means that OCAS becomes equivalent to CPA */
+ /* see paper Franc&Sonneburg JMLR 2009 */
static const uint32_t QPSolverMaxIter = 10000000;
@@ -49,6 +50,477 @@ static double get_time()
return 0.0;
}
+
+/*----------------------------------------------------------------------
+ Linear binary Ocas-SVM solver with additinal contraint enforceing
+ a subset of weights (indices of the weights given by num_nnw/nnw_idx)
+ to be non-negative.
+
+ ----------------------------------------------------------------------*/
+ocas_return_value_T svm_ocas_solver_nnw(
+ double C,
+ uint32_t nData,
+ uint32_t num_nnw,
+ uint32_t* nnw_idx,
+ double TolRel,
+ double TolAbs,
+ double QPBound,
+ double MaxTime,
+ uint32_t _BufSize,
+ uint8_t Method,
+ int (*add_pw_constr)(uint32_t, uint32_t, void*),
+ void (*clip_neg_W)(uint32_t, uint32_t*, void*),
+ void (*compute_W)(double*, double*, double*, uint32_t, void*),
+ double (*update_W)(double, void*),
+ int (*add_new_cut)(double*, uint32_t*, uint32_t, uint32_t, void*),
+ int (*compute_output)(double*, void* ),
+ int (*sort)(double*, double*, uint32_t),
+ void (*ocas_print)(ocas_return_value_T),
+ void* user_data)
+{
+ ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
+ double *b, *alpha, *diag_H;
+ double *output, *old_output;
+ double xi, sq_norm_W, QPSolverTolRel, dot_prod_WoldW, sq_norm_oldW;
+ double A0, B0, GradVal, t, t1, t2, *Ci, *Bi, *hpf, *hpb;
+ double start_time, ocas_start_time;
+ uint32_t cut_length;
+ uint32_t i, *new_cut;
+ uint32_t *I;
+/* uint8_t S = 1; */
+ libqp_state_T qp_exitflag;
+
+ double max_cp_norm;
+ double max_b;
+ double _C[2];
+ uint8_t S[2];
+
+ ocas_start_time = get_time();
+ ocas.qp_solver_time = 0;
+ ocas.output_time = 0;
+ ocas.sort_time = 0;
+ ocas.add_time = 0;
+ ocas.w_time = 0;
+ ocas.print_time = 0;
+
+ BufSize = _BufSize;
+
+ QPSolverTolRel = TolRel*0.5;
+
+ H=NULL;
+ b=NULL;
+ alpha=NULL;
+ new_cut=NULL;
+ I=NULL;
+ diag_H=NULL;
+ output=NULL;
+ old_output=NULL;
+ hpf=NULL;
+ hpb = NULL;
+ Ci=NULL;
+ Bi=NULL;
+
+ /* Hessian matrix contains dot product of normal vectors of selected cutting planes */
+ H = (double*)LIBOCAS_CALLOC(BufSize*BufSize,sizeof(double));
+ if(H == NULL)
+ {
+ ocas.exitflag=-2;
+ goto cleanup;
+ }
+
+ /* bias of cutting planes */
+ b = (double*)LIBOCAS_CALLOC(BufSize,sizeof(double));
+ if(b == NULL)
+ {
+ ocas.exitflag=-2;
+ goto cleanup;
+ }
+
+ alpha = (double*)LIBOCAS_CALLOC(BufSize,sizeof(double));
+ if(alpha == NULL)
+ {
+ ocas.exitflag=-2;
+ goto cleanup;
+ }
+
+ /* indices of examples which define a new cut */
+ new_cut = (uint32_t*)LIBOCAS_CALLOC(nData,sizeof(uint32_t));
+ if(new_cut == NULL)
+ {
+ ocas.exitflag=-2;
+ goto cleanup;
+ }
+
+ I = (uint32_t*)LIBOCAS_CALLOC(BufSize,sizeof(uint32_t));
+ if(I == NULL)
+ {
+ ocas.exitflag=-2;
+ goto cleanup;
+ }
+
+ for(i=0; i< BufSize; i++) I[i] = 2;
+
+ diag_H = (double*)LIBOCAS_CALLOC(BufSize,sizeof(double));
+ if(diag_H == NULL)
+ {
+ ocas.exitflag=-2;
+ goto cleanup;
+ }
+
+ output = (double*)LIBOCAS_CALLOC(nData,sizeof(double));
+ if(output == NULL)
+ {
+ ocas.exitflag=-2;
+ goto cleanup;
+ }
+
+ old_output = (double*)LIBOCAS_CALLOC(nData,sizeof(double));
+ if(old_output == NULL)
+ {
+ ocas.exitflag=-2;
+ goto cleanup;
+ }
+
+ /* array of hinge points used in line-serach */
+ hpf = (double*) LIBOCAS_CALLOC(nData, sizeof(hpf[0]));
+ if(hpf == NULL)
+ {
+ ocas.exitflag=-2;
+ goto cleanup;
+ }
+
+ hpb = (double*) LIBOCAS_CALLOC(nData, sizeof(hpb[0]));
+ if(hpb == NULL)
+ {
+ ocas.exitflag=-2;
+ goto cleanup;
+ }
+
+ /* vectors Ci, Bi are used in the line search procedure */
+ Ci = (double*)LIBOCAS_CALLOC(nData,sizeof(double));
+ if(Ci == NULL)
+ {
+ ocas.exitflag=-2;
+ goto cleanup;
+ }
+
+ Bi = (double*)LIBOCAS_CALLOC(nData,sizeof(double));
+ if(Bi == NULL)
+ {
+ ocas.exitflag=-2;
+ goto cleanup;
+ }
+
+ /* initial cutting planes implementing the non-negativity constraints on W*/
+ _C[0]=10000000.0;
+ _C[1]=C;
+ S[0]=1;
+ S[1]=1;
+ for(i=0; i < num_nnw; i++)
+ {
+ if(add_pw_constr(nnw_idx[i],i, user_data) != 0)
+ {
+ ocas.exitflag=-2;
+ goto cleanup;
+ }
+ diag_H[i] = 1.0;
+ H[LIBOCAS_INDEX(i,i,BufSize)] = 1.0;
+ I[i] = 1;
+ }
+
+ max_cp_norm = 1;
+ max_b = 0;
+
+ /* */
+ ocas.nCutPlanes = num_nnw;
+ ocas.exitflag = 0;
+ ocas.nIter = 0;
+
+ /* Compute initial value of Q_P assuming that W is zero vector.*/
+ sq_norm_W = 0;
+ xi = nData;
+ ocas.Q_P = 0.5*sq_norm_W + C*xi;
+ ocas.Q_D = 0;
+
+ /* Compute the initial cutting plane */
+ cut_length = nData;
+ for(i=0; i < nData; i++)
+ new_cut[i] = i;
+
+ ocas.trn_err = nData;
+ ocas.ocas_time = get_time() - ocas_start_time;
+ /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, Q_P-Q_D/abs(Q_P)=%f\n",
+ ocas.nIter,cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P));
+ */
+ ocas_print(ocas);
+
+ /* main loop */
+ while( ocas.exitflag == 0 )
+ {
+ ocas.nIter++;
+
+ /* append a new cut to the buffer and update H */
+ b[ocas.nCutPlanes] = -(double)cut_length;
+
+ max_b = LIBOCAS_MAX(max_b,(double)cut_length);
+
+ start_time = get_time();
+
+ if(add_new_cut( &H[LIBOCAS_INDEX(0,ocas.nCutPlanes,BufSize)], new_cut, cut_length, ocas.nCutPlanes, user_data ) != 0)
+ {
+ ocas.exitflag=-2;
+ goto cleanup;
+ }
+
+ ocas.add_time += get_time() - start_time;
+
+ /* copy new added row: H(ocas.nCutPlanes,ocas.nCutPlanes,1:ocas.nCutPlanes-1) = H(1:ocas.nCutPlanes-1:ocas.nCutPlanes)' */
+ diag_H[ocas.nCutPlanes] = H[LIBOCAS_INDEX(ocas.nCutPlanes,ocas.nCutPlanes,BufSize)];
+ for(i=0; i < ocas.nCutPlanes; i++) {
+ H[LIBOCAS_INDEX(ocas.nCutPlanes,i,BufSize)] = H[LIBOCAS_INDEX(i,ocas.nCutPlanes,BufSize)];
+ }
+
+ max_cp_norm = LIBOCAS_MAX(max_cp_norm, sqrt(diag_H[ocas.nCutPlanes]));
+
+
+ ocas.nCutPlanes++;
+
+ /* call inner QP solver */
+ start_time = get_time();
+
+ /* compute upper bound on sum of dual variables associated with the positivity constraints */
+ _C[0] = sqrt((double)nData)*(sqrt(C)*sqrt(max_b) + C*max_cp_norm);
+
+/* qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, &C, I, &S, alpha,
+ ocas.nCutPlanes, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBOCAS_PLUS_INF,0);*/
+ qp_exitflag = libqp_splx_solver(&get_col, diag_H, b, _C, I, S, alpha,
+ ocas.nCutPlanes, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBOCAS_PLUS_INF,0);
+
+ ocas.qp_exitflag = qp_exitflag.exitflag;
+
+ ocas.qp_solver_time += get_time() - start_time;
+ ocas.Q_D = -qp_exitflag.QP;
+
+ ocas.nNZAlpha = 0;
+ for(i=0; i < ocas.nCutPlanes; i++) {
+ if( alpha[i] != 0) ocas.nNZAlpha++;
+ }
+
+ sq_norm_oldW = sq_norm_W;
+ start_time = get_time();
+ compute_W( &sq_norm_W, &dot_prod_WoldW, alpha, ocas.nCutPlanes, user_data );
+ clip_neg_W(num_nnw, nnw_idx, user_data);
+ ocas.w_time += get_time() - start_time;
+
+ /* select a new cut */
+ switch( Method )
+ {
+ /* cutting plane algorithm implemented in SVMperf and BMRM */
+ case 0:
+
+ start_time = get_time();
+ if( compute_output( output, user_data ) != 0)
+ {
+ ocas.exitflag=-2;
+ goto cleanup;
+ }
+ ocas.output_time += get_time()-start_time;
+
+ xi = 0;
+ cut_length = 0;
+ ocas.trn_err = 0;
+ for(i=0; i < nData; i++)
+ {
+ if(output[i] <= 0) ocas.trn_err++;
+
+ if(output[i] <= 1) {
+ xi += 1 - output[i];
+ new_cut[cut_length] = i;
+ cut_length++;
+ }
+ }
+ ocas.Q_P = 0.5*sq_norm_W + C*xi;
+
+ ocas.ocas_time = get_time() - ocas_start_time;
+
+ /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, 1-Q_D/Q_P=%f, nza=%4d, err=%.2f%%, qpf=%d\n",
+ ocas.nIter,cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P),
+ ocas.nNZAlpha, 100*(double)ocas.trn_err/(double)nData, ocas.qp_exitflag );
+ */
+
+ start_time = get_time();
+ ocas_print(ocas);
+ ocas.print_time += get_time() - start_time;
+
+ break;
+
+
+ /* Ocas strategy */
+ case 1:
+
+ /* Linesearch */
+ A0 = sq_norm_W -2*dot_prod_WoldW + sq_norm_oldW;
+ B0 = dot_prod_WoldW - sq_norm_oldW;
+
+ memcpy( old_output, output, sizeof(double)*nData );
+
+ start_time = get_time();
+ if( compute_output( output, user_data ) != 0)
+ {
+ ocas.exitflag=-2;
+ goto cleanup;
+ }
+ ocas.output_time += get_time()-start_time;
+
+ uint32_t num_hp = 0;
+ GradVal = B0;
+ for(i=0; i< nData; i++) {
+
+ Ci[i] = C*(1-old_output[i]);
+ Bi[i] = C*(old_output[i] - output[i]);
+
+ double val;
+ if(Bi[i] != 0)
+ val = -Ci[i]/Bi[i];
+ else
+ val = -LIBOCAS_PLUS_INF;
+
+ if (val>0)
+ {
+/* hpi[num_hp] = i;*/
+ hpb[num_hp] = Bi[i];
+ hpf[num_hp] = val;
+ num_hp++;
+ }
+
+ if( (Bi[i] < 0 && val > 0) || (Bi[i] > 0 && val <= 0))
+ GradVal += Bi[i];
+
+ }
+
+ t = 0;
+ if( GradVal < 0 )
+ {
+ start_time = get_time();
+/* if( sort(hpf, hpi, num_hp) != 0)*/
+ if( sort(hpf, hpb, num_hp) != 0 )
+ {
+ ocas.exitflag=-2;
+ goto cleanup;
+ }
+ ocas.sort_time += get_time() - start_time;
+
+ double t_new, GradVal_new;
+ i = 0;
+ while( GradVal < 0 && i < num_hp )
+ {
+ t_new = hpf[i];
+ GradVal_new = GradVal + LIBOCAS_ABS(hpb[i]) + A0*(t_new-t);
+
+ if( GradVal_new >= 0 )
+ {
+ t = t + GradVal*(t-t_new)/(GradVal_new - GradVal);
+ }
+ else
+ {
+ t = t_new;
+ i++;
+ }
+
+ GradVal = GradVal_new;
+ }
+ }
+
+ /*
+ t = hpf[0] - 1;
+ i = 0;
+ GradVal = t*A0 + Bsum;
+ while( GradVal < 0 && i < num_hp && hpf[i] < LIBOCAS_PLUS_INF ) {
+ t = hpf[i];
+ Bsum = Bsum + LIBOCAS_ABS(Bi[hpi[i]]);
+ GradVal = t*A0 + Bsum;
+ i++;
+ }
+ */
+ t = LIBOCAS_MAX(t,0); /* just sanity check; t < 0 should not ocure */
+
+ /* this guarantees that the new solution will not violate the positivity constraints on W */
+ t = LIBOCAS_MIN(t,1);
+
+ t1 = t; /* new (best so far) W */
+ t2 = t+MU*(1.0-t); /* new cutting plane */
+ /* t2 = t+(1.0-t)/10.0; */
+
+ /* update W to be the best so far solution */
+ sq_norm_W = update_W( t1, user_data );
+
+ /* select a new cut */
+ xi = 0;
+ cut_length = 0;
+ ocas.trn_err = 0;
+ for(i=0; i < nData; i++ ) {
+
+ if( (old_output[i]*(1-t2) + t2*output[i]) <= 1 )
+ {
+ new_cut[cut_length] = i;
+ cut_length++;
+ }
+
+ output[i] = old_output[i]*(1-t1) + t1*output[i];
+
+ if( output[i] <= 1) xi += 1-output[i];
+ if( output[i] <= 0) ocas.trn_err++;
+
+ }
+
+ ocas.Q_P = 0.5*sq_norm_W + C*xi;
+
+ ocas.ocas_time = get_time() - ocas_start_time;
+
+ /* ocas_print("%4d: tim=%f, Q_P=%f, Q_D=%f, Q_P-Q_D=%f, 1-Q_D/Q_P=%f, nza=%4d, err=%.2f%%, qpf=%d\n",
+ ocas.nIter, cur_time, ocas.Q_P,ocas.Q_D,ocas.Q_P-ocas.Q_D,(ocas.Q_P-ocas.Q_D)/LIBOCAS_ABS(ocas.Q_P),
+ ocas.nNZAlpha, 100*(double)ocas.trn_err/(double)nData, ocas.qp_exitflag );
+ */
+
+ start_time = get_time();
+ ocas_print(ocas);
+ ocas.print_time += get_time() - start_time;
+
+ break;
+ }
+
+ /* Stopping conditions */
+ if( ocas.Q_P - ocas.Q_D <= TolRel*LIBOCAS_ABS(ocas.Q_P)) ocas.exitflag = 1;
+ if( ocas.Q_P - ocas.Q_D <= TolAbs) ocas.exitflag = 2;
+ if( ocas.Q_P <= QPBound) ocas.exitflag = 3;
+ if( ocas.ocas_time >= MaxTime) ocas.exitflag = 4;
+ if(ocas.nCutPlanes >= BufSize) ocas.exitflag = -1;
+
+ } /* end of the main loop */
+
+cleanup:
+
+ LIBOCAS_FREE(H);
+ LIBOCAS_FREE(b);
+ LIBOCAS_FREE(alpha);
+ LIBOCAS_FREE(new_cut);
+ LIBOCAS_FREE(I);
+ LIBOCAS_FREE(diag_H);
+ LIBOCAS_FREE(output);
+ LIBOCAS_FREE(old_output);
+ LIBOCAS_FREE(hpf);
+/* LIBOCAS_FREE(hpi);*/
+ LIBOCAS_FREE(hpb);
+ LIBOCAS_FREE(Ci);
+ LIBOCAS_FREE(Bi);
+
+ ocas.ocas_time = get_time() - ocas_start_time;
+
+ return(ocas);
+}
+
+
+
/*----------------------------------------------------------------------
Linear binary Ocas-SVM solver.
----------------------------------------------------------------------*/
@@ -69,7 +541,7 @@ ocas_return_value_T svm_ocas_solver(
void (*ocas_print)(ocas_return_value_T),
void* user_data)
{
- ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
+ ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
double *b, *alpha, *diag_H;
double *output, *old_output;
double xi, sq_norm_W, QPSolverTolRel, dot_prod_WoldW, sq_norm_oldW;
@@ -399,7 +871,7 @@ ocas_return_value_T svm_ocas_solver(
t = LIBOCAS_MAX(t,0); /* just sanity check; t < 0 should not ocure */
t1 = t; /* new (best so far) W */
- t2 = t+LAMBDA*(1.0-t); /* new cutting plane */
+ t2 = t+MU*(1.0-t); /* new cutting plane */
/* t2 = t+(1.0-t)/10.0; */
/* update W to be the best so far solution */
@@ -492,7 +964,7 @@ ocas_return_value_T svm_ocas_solver_difC(
void (*ocas_print)(ocas_return_value_T),
void* user_data)
{
- ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
+ ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
double *b, *alpha, *diag_H;
double *output, *old_output;
double xi, sq_norm_W, QPSolverTolRel, dot_prod_WoldW, sq_norm_oldW;
@@ -841,7 +1313,7 @@ ocas_return_value_T svm_ocas_solver_difC(
t = LIBOCAS_MAX(t,0); /* just sanity check; t < 0 should not ocure */
t1 = t; /* new (best so far) W */
- t2 = t+(1.0-t)*LAMBDA; /* new cutting plane */
+ t2 = t+(1.0-t)*MU; /* new cutting plane */
/* t2 = t+(1.0-t)/10.0; new cutting plane */
/* update W to be the best so far solution */
@@ -938,7 +1410,7 @@ static void findactive(double *Theta, double *SortedA, uint32_t *nSortedA, doubl
tmp = B[0];
idx = 0;
i = 0;
- while( i < n-1 && A[i] == A[i+1])
+ while( i < (uint32_t)n-1 && A[i] == A[i+1])
{
if( B[i+1] > B[idx] )
{
@@ -954,11 +1426,11 @@ static void findactive(double *Theta, double *SortedA, uint32_t *nSortedA, doubl
while(1)
{
start = idx + 1;
- while( start < n && A[idx] == A[start])
+ while( start < (uint32_t)n && A[idx] == A[start])
start++;
theta = LIBOCAS_PLUS_INF;
- for(j=start; j < n; j++)
+ for(j=start; j < (uint32_t)n; j++)
{
tmp = (B[j] - B[idx])/(A[idx]-A[j]);
if( tmp < theta)
@@ -1003,7 +1475,7 @@ ocas_return_value_T msvm_ocas_solver(
void (*ocas_print)(ocas_return_value_T),
void* user_data)
{
- ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
+ ocas_return_value_T ocas={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
double *b, *alpha, *diag_H;
double *output, *old_output;
double xi, sq_norm_W, QPSolverTolRel, QPSolverTolAbs, dot_prod_WoldW, sq_norm_oldW;
@@ -1384,7 +1856,7 @@ ocas_return_value_T msvm_ocas_solver(
t = min_x;
t1 = t; /* new (best so far) W */
- t2 = t+(1.0-t)*LAMBDA; /* new cutting plane */
+ t2 = t+(1.0-t)*MU; /* new cutting plane */
/* t2 = t+(1.0-t)/10.0; */
/* update W to be the best so far solution */
diff --git a/libocas.h b/libocas.h
index f1b45ed..1ea1ed8 100644
--- a/libocas.h
+++ b/libocas.h
@@ -119,5 +119,29 @@ ocas_return_value_T msvm_ocas_solver(
void (*ocas_print)(ocas_return_value_T),
void* user_data);
+
+/* binary linear SVM solver */
+ocas_return_value_T svm_ocas_solver_nnw(
+ double C, /* regularizarion constant */
+ uint32_t nData, /* number of exmaples */
+ uint32_t num_nnw, /* number of components of W which must non-negative*/
+ uint32_t* nnw_idx, /* indices of W which must be non-negative */
+ double TolRel, /* halts if 1-Q_P/Q_D <= TolRel */
+ double TolAbs, /* halts if Q_P-Q_D <= TolRel */
+ double QPBound, /* halts if QP <= QPBound */
+ double MaxTime, /* maximal time in seconds spent in optmization */
+ uint32_t BufSize, /* maximal number of buffered cutting planes */
+ uint8_t Method, /* 0..standard CP (SVM-Perf,BMRM), 1..OCAS */
+ int (*add_pw_constr)(uint32_t, uint32_t, void*),
+ void (*clip_neg_w)(uint32_t, uint32_t*, void*),
+ void (*compute_W)(double*, double*, double*, uint32_t, void*),
+ double (*update_W)(double, void*),
+ int (*add_new_cut)(double*, uint32_t*, uint32_t, uint32_t, void*),
+ int (*compute_output)( double*, void* ),
+ int (*sort)(double*, double*, uint32_t),
+ void (*ocas_print)(ocas_return_value_T),
+ void* user_data);
+
+
#endif /* libocas_h */
diff --git a/libocas.so b/libocas.so
deleted file mode 100755
index 41a9953..0000000
Binary files a/libocas.so and /dev/null differ
diff --git a/libqp_splx.c b/libqp_splx.c
index 37a36eb..3c673e3 100644
--- a/libqp_splx.c
+++ b/libqp_splx.c
@@ -1,407 +1,407 @@
-/*-----------------------------------------------------------------------
- * libqp_splx.c: solver for Quadratic Programming task with
- * simplex constraints.
- *
- * DESCRIPTION
- * The library provides function which solves the following instance of
- * a convex Quadratic Programmin task:
- *
- * min QP(x):= 0.5*x'*H*x + f'*x
- * x
- *
- * subject to:
- * sum_{i in I_k} x[i] == b[k] for all k such that S[k] == 0
- * sum_{i in I_k} x[i] <= b[k] for all k such that S[k] == 1
- * x(i) >= 0 for all i=1:n
- *
- * where I_k = { i | I[i] == k}, k={1,...,m}.
- *
- * A precision of the found solution is controled by the input argumens
- * MaxIter, TolAbs, QP_TH and MaxIter which define the stopping conditions:
- *
- * nIter >= MaxIter -> exitflag = 0 Number of iterations
- * QP-QD <= TolAbs -> exitflag = 1 Abs. tolerance (duality gap)
- * QP-QD <= QP*TolRel -> exitflag = 2 Relative tolerance
- * QP <= QP_TH -> exitflag = 3 Threshold on objective value
- *
- * where QP and QD are primal respectively dual objective values.
- *
- * INPUT ARGUMENTS
- * get_col function which returns pointer to the i-th column of H.
- * diag_H [double n x 1] vector containing values on the diagonal of H.
- * f [double n x 1] vector.
- * b [double n x 1] vector of positive numbers.
- * I [uint16_T n x 1] vector containing numbers 1...m.
- * S [uint8_T n x 1] vector containing numbers 0 and 1.
- * x [double n x 1] solution vector; must be feasible.
- * n [uint32_t 1 x 1] dimension of H.
- * MaxIter [uint32_t 1 x 1] max number of iterations.
- * TolAbs [double 1 x 1] Absolute tolerance.
- * TolRel [double 1 x 1] Relative tolerance.
- * QP_TH [double 1 x 1] Threshold on the primal value.
- * print_state print function; if == NULL it is not called.
- *
- * RETURN VALUE
- * structure [libqp_state_T]
- * .QP [1 x 1] Primal objective value.
- * .QD [1 x 1] Dual objective value.
- * .nIter [1 x 1] Number of iterations.
- * .exitflag [1 x 1] Indicates which stopping condition was used:
- * -1 ... Not enough memory.
- * 0 ... Maximal number of iteations reached: nIter >= MaxIter.
- * 1 ... Relarive tolerance reached: QP-QD <= abs(QP)*TolRel
- * 2 ... Absolute tolerance reached: QP-QD <= TolAbs
- * 3 ... Objective value reached threshold: QP <= QP_TH.
- *
- * REFERENCE
- * The algorithm is described in:
- * V. Franc, V. Hlavac. A Novel Algorithm for Learning Support Vector Machines
- * with Structured Output Spaces. Research Report K333 22/06, CTU-CMP-2006-04.
- * May, 2006. ftp://cmp.felk.cvut.cz/pub/cmp/articles/franc/Franc-TR-2006-04.ps
- *
- * Copyright (C) 2006-2008 Vojtech Franc, xfrancv at cmp.felk.cvut.cz
- * Center for Machine Perception, CTU FEL Prague
- *
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public
- * License as published by the Free Software Foundation;
- * Version 3, 29 June 2007
- *-------------------------------------------------------------------- */
-
-#include <math.h>
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-#include <stdint.h>
-#include <limits.h>
-
-#include "libqp.h"
-
-libqp_state_T libqp_splx_solver(const double* (*get_col)(uint32_t),
- double *diag_H,
- double *f,
- double *b,
- uint32_t *I,
- uint8_t *S,
- double *x,
- uint32_t n,
- uint32_t MaxIter,
- double TolAbs,
- double TolRel,
- double QP_TH,
- void (*print_state)(libqp_state_T state))
-{
- double *d;
- double *col_u, *col_v;
- double *x_neq;
- double tmp;
- double improv;
- double tmp_num;
- double tmp_den=0;
- double tau=0;
- double delta;
- uint32_t *inx;
- uint32_t *nk;
- uint32_t m;
- uint32_t u=0;
- uint32_t v=0;
- uint32_t k;
- uint32_t i, j;
- libqp_state_T state;
-
-
- /* ------------------------------------------------------------
- Initialization
- ------------------------------------------------------------ */
- state.nIter = 0;
- state.QP = LIBQP_PLUS_INF;
- state.QD = -LIBQP_PLUS_INF;
- state.exitflag = 100;
-
- inx=NULL;
- nk=NULL;
- d=NULL;
- x_neq = NULL;
-
- /* count number of constraints */
- for( i=0, m=0; i < n; i++ )
- m = LIBQP_MAX(m,I[i]);
-
- /* auxciliary variables for tranforming equalities to inequalities */
- x_neq = (double*) LIBQP_CALLOC(m, sizeof(double));
- if( x_neq == NULL )
- {
- state.exitflag=-1;
- goto cleanup;
- }
-
- /* inx is translation table between variable index i and its contraint */
- inx = (uint32_t*) LIBQP_CALLOC(m*n, sizeof(uint32_t));
- if( inx == NULL )
- {
- state.exitflag=-1;
- goto cleanup;
- }
-
- /* nk is the number of variables coupled by i-th linear constraint */
- nk = (uint32_t*) LIBQP_CALLOC(m, sizeof(uint32_t));
- if( nk == NULL )
- {
- state.exitflag=-1;
- goto cleanup;
- }
-
- /* setup auxciliary variables */
- for( i=0; i < m; i++ )
- x_neq[i] = b[i];
-
-
- /* create inx and nk */
- for( i=0; i < n; i++ ) {
- k = I[i]-1;
- inx[LIBQP_INDEX(nk[k],k,n)] = i;
- nk[k]++;
-
- if(S[k] != 0)
- x_neq[k] -= x[i];
- }
-
- /* d = H*x + f is gradient*/
- d = (double*) LIBQP_CALLOC(n, sizeof(double));
- if( d == NULL )
- {
- state.exitflag=-1;
- goto cleanup;
- }
-
- /* compute gradient */
- for( i=0; i < n; i++ )
- {
- d[i] += f[i];
- if( x[i] > 0 ) {
- col_u = (double*)get_col(i);
- for( j=0; j < n; j++ ) {
- d[j] += col_u[j]*x[i];
- }
- }
- }
-
- /* compute state.QP = 0.5*x'*(f+d);
- state.QD = 0.5*x'*(f-d); */
- for( i=0, state.QP = 0, state.QD=0; i < n; i++)
- {
- state.QP += x[i]*(f[i]+d[i]);
- state.QD += x[i]*(f[i]-d[i]);
- }
- state.QP = 0.5*state.QP;
- state.QD = 0.5*state.QD;
-
- for( i=0; i < m; i++ )
- {
- for( j=0, tmp = LIBQP_PLUS_INF; j < nk[i]; j++ )
- tmp = LIBQP_MIN(tmp, d[inx[LIBQP_INDEX(j,i,n)]]);
-
- if(S[i] == 0)
- state.QD += b[i]*tmp;
- else
- state.QD += b[i]*LIBQP_MIN(tmp,0);
- }
-
- /* print initial state */
- if( print_state != NULL)
- print_state( state );
-
- /* ------------------------------------------------------------
- Main optimization loop
- ------------------------------------------------------------ */
- while( state.exitflag == 100 )
- {
- state.nIter ++;
-
- /* go over blocks of variables coupled by lin. constraint */
- for( k=0; k < m; k++ )
- {
-
- /* compute u = argmin_{i in I_k} d[i]
- delta = sum_{i in I_k} x[i]*d[i] - b*min_{i in I_k} */
- for( j=0, tmp = LIBQP_PLUS_INF, delta = 0; j < nk[k]; j++ )
- {
- i = inx[LIBQP_INDEX(j,k,n)];
- delta += x[i]*d[i];
- if( tmp > d[i] ) {
- tmp = d[i];
- u = i;
- }
- }
-
- if(S[k] != 0 && d[u] > 0)
- u = -1;
- else
- delta -= b[k]*d[u];
-
- /* if satisfied then k-th block of variables needs update */
- if( delta > TolAbs/m && delta > TolRel*LIBQP_ABS(state.QP)/m)
- {
- /* for fixed u select v = argmax_{i in I_k} Improvement(i) */
- if( u != -1 )
- {
- col_u = (double*)get_col(u);
- improv = -LIBQP_PLUS_INF;
- for( j=0; j < nk[k]; j++ )
- {
- i = inx[LIBQP_INDEX(j,k,n)];
-
- if(x[i] > 0 && i != u)
- {
- tmp_num = x[i]*(d[i] - d[u]);
- tmp_den = x[i]*x[i]*(diag_H[u] - 2*col_u[i] + diag_H[i]);
- if( tmp_den > 0 )
- {
- if( tmp_num < tmp_den )
- tmp = tmp_num*tmp_num / tmp_den;
- else
- tmp = tmp_num - 0.5 * tmp_den;
-
- if( tmp > improv )
- {
- improv = tmp;
- tau = LIBQP_MIN(1,tmp_num/tmp_den);
- v = i;
- }
- }
- }
- }
-
- /* check if virtual variable can be for updated */
- if(x_neq[k] > 0 && S[k] != 0)
- {
- tmp_num = -x_neq[k]*d[u];
- tmp_den = x_neq[k]*x_neq[k]*diag_H[u];
- if( tmp_den > 0 )
- {
- if( tmp_num < tmp_den )
- tmp = tmp_num*tmp_num / tmp_den;
- else
- tmp = tmp_num - 0.5 * tmp_den;
-
- if( tmp > improv )
- {
- improv = tmp;
- tau = LIBQP_MIN(1,tmp_num/tmp_den);
- v = -1;
- }
- }
- }
-
- /* minimize objective w.r.t variable u and v */
- if(v != -1)
- {
- tmp = x[v]*tau;
- x[u] += tmp;
- x[v] -= tmp;
-
- /* update d = H*x + f */
- col_v = (double*)get_col(v);
- for(i = 0; i < n; i++ )
- d[i] += tmp*(col_u[i]-col_v[i]);
- }
- else
- {
- tmp = x_neq[k]*tau;
- x[u] += tmp;
- x_neq[k] -= tmp;
-
- /* update d = H*x + f */
- for(i = 0; i < n; i++ )
- d[i] += tmp*col_u[i];
- }
- }
- else
- {
- improv = -LIBQP_PLUS_INF;
- for( j=0; j < nk[k]; j++ )
- {
- i = inx[LIBQP_INDEX(j,k,n)];
-
- if(x[i] > 0)
- {
- tmp_num = x[i]*d[i];
- tmp_den = x[i]*x[i]*diag_H[i];
- if( tmp_den > 0 )
- {
- if( tmp_num < tmp_den )
- tmp = tmp_num*tmp_num / tmp_den;
- else
- tmp = tmp_num - 0.5 * tmp_den;
-
- if( tmp > improv )
- {
- improv = tmp;
- tau = LIBQP_MIN(1,tmp_num/tmp_den);
- v = i;
- }
- }
- }
- }
-
- tmp = x[v]*tau;
- x_neq[k] += tmp;
- x[v] -= tmp;
-
- /* update d = H*x + f */
- col_v = (double*)get_col(v);
- for(i = 0; i < n; i++ )
- d[i] -= tmp*col_v[i];
- }
-
- /* update objective value */
- state.QP = state.QP - improv;
- }
- }
-
- /* Compute primal and dual objectives */
- for( i=0, state.QP = 0, state.QD=0; i < n; i++)
- {
- state.QP += x[i]*(f[i]+d[i]);
- state.QD += x[i]*(f[i]-d[i]);
- }
- state.QP = 0.5*state.QP;
- state.QD = 0.5*state.QD;
-
- for( k=0; k < m; k++ )
- {
- for( j=0,tmp = LIBQP_PLUS_INF; j < nk[k]; j++ ) {
- i = inx[LIBQP_INDEX(j,k,n)];
- tmp = LIBQP_MIN(tmp, d[i]);
- }
-
- if(S[k] == 0)
- state.QD += b[k]*tmp;
- else
- state.QD += b[k]*LIBQP_MIN(tmp,0);
- }
-
- /* print state */
- if( print_state != NULL)
- print_state( state );
-
- /* check stopping conditions */
- if(state.QP-state.QD <= LIBQP_ABS(state.QP)*TolRel ) state.exitflag = 1;
- else if( state.QP-state.QD <= TolAbs ) state.exitflag = 2;
- else if( state.QP <= QP_TH ) state.exitflag = 3;
- else if( state.nIter >= MaxIter) state.exitflag = 0;
- }
-
- /*----------------------------------------------------------
- Clean up
- ---------------------------------------------------------- */
-cleanup:
- LIBQP_FREE( d );
- LIBQP_FREE( inx );
- LIBQP_FREE( nk );
- LIBQP_FREE( x_neq );
-
- return( state );
-}
-
-
+/*-----------------------------------------------------------------------
+ * libqp_splx.c: solver for Quadratic Programming task with
+ * simplex constraints.
+ *
+ * DESCRIPTION
+ * The library provides function which solves the following instance of
+ * a convex Quadratic Programmin task:
+ *
+ * min QP(x):= 0.5*x'*H*x + f'*x
+ * x
+ *
+ * subject to:
+ * sum_{i in I_k} x[i] == b[k] for all k such that S[k] == 0
+ * sum_{i in I_k} x[i] <= b[k] for all k such that S[k] == 1
+ * x(i) >= 0 for all i=1:n
+ *
+ * where I_k = { i | I[i] == k}, k={1,...,m}.
+ *
+ * A precision of the found solution is controled by the input argumens
+ * MaxIter, TolAbs, QP_TH and MaxIter which define the stopping conditions:
+ *
+ * nIter >= MaxIter -> exitflag = 0 Number of iterations
+ * QP-QD <= TolAbs -> exitflag = 1 Abs. tolerance (duality gap)
+ * QP-QD <= QP*TolRel -> exitflag = 2 Relative tolerance
+ * QP <= QP_TH -> exitflag = 3 Threshold on objective value
+ *
+ * where QP and QD are primal respectively dual objective values.
+ *
+ * INPUT ARGUMENTS
+ * get_col function which returns pointer to the i-th column of H.
+ * diag_H [double n x 1] vector containing values on the diagonal of H.
+ * f [double n x 1] vector.
+ * b [double n x 1] vector of positive numbers.
+ * I [uint16_T n x 1] vector containing numbers 1...m.
+ * S [uint8_T n x 1] vector containing numbers 0 and 1.
+ * x [double n x 1] solution vector; must be feasible.
+ * n [uint32_t 1 x 1] dimension of H.
+ * MaxIter [uint32_t 1 x 1] max number of iterations.
+ * TolAbs [double 1 x 1] Absolute tolerance.
+ * TolRel [double 1 x 1] Relative tolerance.
+ * QP_TH [double 1 x 1] Threshold on the primal value.
+ * print_state print function; if == NULL it is not called.
+ *
+ * RETURN VALUE
+ * structure [libqp_state_T]
+ * .QP [1 x 1] Primal objective value.
+ * .QD [1 x 1] Dual objective value.
+ * .nIter [1 x 1] Number of iterations.
+ * .exitflag [1 x 1] Indicates which stopping condition was used:
+ * -1 ... Not enough memory.
+ * 0 ... Maximal number of iteations reached: nIter >= MaxIter.
+ * 1 ... Relarive tolerance reached: QP-QD <= abs(QP)*TolRel
+ * 2 ... Absolute tolerance reached: QP-QD <= TolAbs
+ * 3 ... Objective value reached threshold: QP <= QP_TH.
+ *
+ * REFERENCE
+ * The algorithm is described in:
+ * V. Franc, V. Hlavac. A Novel Algorithm for Learning Support Vector Machines
+ * with Structured Output Spaces. Research Report K333 22/06, CTU-CMP-2006-04.
+ * May, 2006. ftp://cmp.felk.cvut.cz/pub/cmp/articles/franc/Franc-TR-2006-04.ps
+ *
+ * Copyright (C) 2006-2008 Vojtech Franc, xfrancv at cmp.felk.cvut.cz
+ * Center for Machine Perception, CTU FEL Prague
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public
+ * License as published by the Free Software Foundation;
+ * Version 3, 29 June 2007
+ *-------------------------------------------------------------------- */
+
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <stdint.h>
+#include <limits.h>
+
+#include "libqp.h"
+
+libqp_state_T libqp_splx_solver(const double* (*get_col)(uint32_t),
+ double *diag_H,
+ double *f,
+ double *b,
+ uint32_t *I,
+ uint8_t *S,
+ double *x,
+ uint32_t n,
+ uint32_t MaxIter,
+ double TolAbs,
+ double TolRel,
+ double QP_TH,
+ void (*print_state)(libqp_state_T state))
+{
+ double *d;
+ double *col_u, *col_v;
+ double *x_neq;
+ double tmp;
+ double improv;
+ double tmp_num;
+ double tmp_den=0;
+ double tau=0;
+ double delta;
+ uint32_t *inx;
+ uint32_t *nk;
+ uint32_t m;
+ uint32_t u=0;
+ uint32_t v=0;
+ uint32_t k;
+ uint32_t i, j;
+ libqp_state_T state;
+
+
+ /* ------------------------------------------------------------
+ Initialization
+ ------------------------------------------------------------ */
+ state.nIter = 0;
+ state.QP = LIBQP_PLUS_INF;
+ state.QD = -LIBQP_PLUS_INF;
+ state.exitflag = 100;
+
+ inx=NULL;
+ nk=NULL;
+ d=NULL;
+ x_neq = NULL;
+
+ /* count number of constraints */
+ for( i=0, m=0; i < n; i++ )
+ m = LIBQP_MAX(m,I[i]);
+
+ /* auxciliary variables for tranforming equalities to inequalities */
+ x_neq = (double*) LIBQP_CALLOC(m, sizeof(double));
+ if( x_neq == NULL )
+ {
+ state.exitflag=-1;
+ goto cleanup;
+ }
+
+ /* inx is translation table between variable index i and its contraint */
+ inx = (uint32_t*) LIBQP_CALLOC(m*n, sizeof(uint32_t));
+ if( inx == NULL )
+ {
+ state.exitflag=-1;
+ goto cleanup;
+ }
+
+ /* nk is the number of variables coupled by i-th linear constraint */
+ nk = (uint32_t*) LIBQP_CALLOC(m, sizeof(uint32_t));
+ if( nk == NULL )
+ {
+ state.exitflag=-1;
+ goto cleanup;
+ }
+
+ /* setup auxciliary variables */
+ for( i=0; i < m; i++ )
+ x_neq[i] = b[i];
+
+
+ /* create inx and nk */
+ for( i=0; i < n; i++ ) {
+ k = I[i]-1;
+ inx[LIBQP_INDEX(nk[k],k,n)] = i;
+ nk[k]++;
+
+ if(S[k] != 0)
+ x_neq[k] -= x[i];
+ }
+
+ /* d = H*x + f is gradient*/
+ d = (double*) LIBQP_CALLOC(n, sizeof(double));
+ if( d == NULL )
+ {
+ state.exitflag=-1;
+ goto cleanup;
+ }
+
+ /* compute gradient */
+ for( i=0; i < n; i++ )
+ {
+ d[i] += f[i];
+ if( x[i] > 0 ) {
+ col_u = (double*)get_col(i);
+ for( j=0; j < n; j++ ) {
+ d[j] += col_u[j]*x[i];
+ }
+ }
+ }
+
+ /* compute state.QP = 0.5*x'*(f+d);
+ state.QD = 0.5*x'*(f-d); */
+ for( i=0, state.QP = 0, state.QD=0; i < n; i++)
+ {
+ state.QP += x[i]*(f[i]+d[i]);
+ state.QD += x[i]*(f[i]-d[i]);
+ }
+ state.QP = 0.5*state.QP;
+ state.QD = 0.5*state.QD;
+
+ for( i=0; i < m; i++ )
+ {
+ for( j=0, tmp = LIBQP_PLUS_INF; j < nk[i]; j++ )
+ tmp = LIBQP_MIN(tmp, d[inx[LIBQP_INDEX(j,i,n)]]);
+
+ if(S[i] == 0)
+ state.QD += b[i]*tmp;
+ else
+ state.QD += b[i]*LIBQP_MIN(tmp,0);
+ }
+
+ /* print initial state */
+ if( print_state != NULL)
+ print_state( state );
+
+ /* ------------------------------------------------------------
+ Main optimization loop
+ ------------------------------------------------------------ */
+ while( state.exitflag == 100 )
+ {
+ state.nIter ++;
+
+ /* go over blocks of variables coupled by lin. constraint */
+ for( k=0; k < m; k++ )
+ {
+
+ /* compute u = argmin_{i in I_k} d[i]
+ delta = sum_{i in I_k} x[i]*d[i] - b*min_{i in I_k} */
+ for( j=0, tmp = LIBQP_PLUS_INF, delta = 0; j < nk[k]; j++ )
+ {
+ i = inx[LIBQP_INDEX(j,k,n)];
+ delta += x[i]*d[i];
+ if( tmp > d[i] ) {
+ tmp = d[i];
+ u = i;
+ }
+ }
+
+ if(S[k] != 0 && d[u] > 0)
+ u = -1;
+ else
+ delta -= b[k]*d[u];
+
+ /* if satisfied then k-th block of variables needs update */
+ if( delta > TolAbs/m && delta > TolRel*LIBQP_ABS(state.QP)/m)
+ {
+ /* for fixed u select v = argmax_{i in I_k} Improvement(i) */
+ if( u != -1 )
+ {
+ col_u = (double*)get_col(u);
+ improv = -LIBQP_PLUS_INF;
+ for( j=0; j < nk[k]; j++ )
+ {
+ i = inx[LIBQP_INDEX(j,k,n)];
+
+ if(x[i] > 0 && i != u)
+ {
+ tmp_num = x[i]*(d[i] - d[u]);
+ tmp_den = x[i]*x[i]*(diag_H[u] - 2*col_u[i] + diag_H[i]);
+ if( tmp_den > 0 )
+ {
+ if( tmp_num < tmp_den )
+ tmp = tmp_num*tmp_num / tmp_den;
+ else
+ tmp = tmp_num - 0.5 * tmp_den;
+
+ if( tmp > improv )
+ {
+ improv = tmp;
+ tau = LIBQP_MIN(1,tmp_num/tmp_den);
+ v = i;
+ }
+ }
+ }
+ }
+
+ /* check if virtual variable can be for updated */
+ if(x_neq[k] > 0 && S[k] != 0)
+ {
+ tmp_num = -x_neq[k]*d[u];
+ tmp_den = x_neq[k]*x_neq[k]*diag_H[u];
+ if( tmp_den > 0 )
+ {
+ if( tmp_num < tmp_den )
+ tmp = tmp_num*tmp_num / tmp_den;
+ else
+ tmp = tmp_num - 0.5 * tmp_den;
+
+ if( tmp > improv )
+ {
+ improv = tmp;
+ tau = LIBQP_MIN(1,tmp_num/tmp_den);
+ v = -1;
+ }
+ }
+ }
+
+ /* minimize objective w.r.t variable u and v */
+ if(v != -1)
+ {
+ tmp = x[v]*tau;
+ x[u] += tmp;
+ x[v] -= tmp;
+
+ /* update d = H*x + f */
+ col_v = (double*)get_col(v);
+ for(i = 0; i < n; i++ )
+ d[i] += tmp*(col_u[i]-col_v[i]);
+ }
+ else
+ {
+ tmp = x_neq[k]*tau;
+ x[u] += tmp;
+ x_neq[k] -= tmp;
+
+ /* update d = H*x + f */
+ for(i = 0; i < n; i++ )
+ d[i] += tmp*col_u[i];
+ }
+ }
+ else
+ {
+ improv = -LIBQP_PLUS_INF;
+ for( j=0; j < nk[k]; j++ )
+ {
+ i = inx[LIBQP_INDEX(j,k,n)];
+
+ if(x[i] > 0)
+ {
+ tmp_num = x[i]*d[i];
+ tmp_den = x[i]*x[i]*diag_H[i];
+ if( tmp_den > 0 )
+ {
+ if( tmp_num < tmp_den )
+ tmp = tmp_num*tmp_num / tmp_den;
+ else
+ tmp = tmp_num - 0.5 * tmp_den;
+
+ if( tmp > improv )
+ {
+ improv = tmp;
+ tau = LIBQP_MIN(1,tmp_num/tmp_den);
+ v = i;
+ }
+ }
+ }
+ }
+
+ tmp = x[v]*tau;
+ x_neq[k] += tmp;
+ x[v] -= tmp;
+
+ /* update d = H*x + f */
+ col_v = (double*)get_col(v);
+ for(i = 0; i < n; i++ )
+ d[i] -= tmp*col_v[i];
+ }
+
+ /* update objective value */
+ state.QP = state.QP - improv;
+ }
+ }
+
+ /* Compute primal and dual objectives */
+ for( i=0, state.QP = 0, state.QD=0; i < n; i++)
+ {
+ state.QP += x[i]*(f[i]+d[i]);
+ state.QD += x[i]*(f[i]-d[i]);
+ }
+ state.QP = 0.5*state.QP;
+ state.QD = 0.5*state.QD;
+
+ for( k=0; k < m; k++ )
+ {
+ for( j=0,tmp = LIBQP_PLUS_INF; j < nk[k]; j++ ) {
+ i = inx[LIBQP_INDEX(j,k,n)];
+ tmp = LIBQP_MIN(tmp, d[i]);
+ }
+
+ if(S[k] == 0)
+ state.QD += b[k]*tmp;
+ else
+ state.QD += b[k]*LIBQP_MIN(tmp,0);
+ }
+
+ /* print state */
+ if( print_state != NULL)
+ print_state( state );
+
+ /* check stopping conditions */
+ if(state.QP-state.QD <= LIBQP_ABS(state.QP)*TolRel ) state.exitflag = 1;
+ else if( state.QP-state.QD <= TolAbs ) state.exitflag = 2;
+ else if( state.QP <= QP_TH ) state.exitflag = 3;
+ else if( state.nIter >= MaxIter) state.exitflag = 0;
+ }
+
+ /*----------------------------------------------------------
+ Clean up
+ ---------------------------------------------------------- */
+cleanup:
+ LIBQP_FREE( d );
+ LIBQP_FREE( inx );
+ LIBQP_FREE( nk );
+ LIBQP_FREE( x_neq );
+
+ return( state );
+}
+
+
diff --git a/linclassif b/linclassif
deleted file mode 100755
index b5c20eb..0000000
Binary files a/linclassif and /dev/null differ
diff --git a/linclassif_light.mexa64 b/linclassif_light.mexa64
deleted file mode 100755
index a13e38f..0000000
Binary files a/linclassif_light.mexa64 and /dev/null differ
diff --git a/load_svmlight_file.c b/load_svmlight_file.c
new file mode 100644
index 0000000..691d9fb
--- /dev/null
+++ b/load_svmlight_file.c
@@ -0,0 +1,244 @@
+/*=================================================================
+ * [feat,labels] = load_svmlight_format(file_name)
+ * [feat,labels] = load_svmlight_format(file_name,verb)
+ *
+ * This function reads examples from a file complaying to SVM^light
+ * format.
+ *
+ *
+ *
+ *
+ *=================================================================*/
+
+#include <stdio.h>
+#include <string.h>
+#include <stdint.h>
+#include <mex.h>
+#include <sys/time.h>
+#include <time.h>
+#include <errno.h>
+
+#if !defined(MX_API_VER) || MX_API_VER<0x07040000
+#define mwSize int
+#define INDEX_TYPE_T int
+#define mwIndex int
+#else
+#define INDEX_TYPE_T mwSize
+#endif
+
+#include "lib_svmlight_format.h"
+
+#define MaxExamples 50000000
+
+#define MIN(A,B) ((A) > (B) ? (B) : (A))
+#define MAX(A,B) ((A) < (B) ? (B) : (A))
+#define ABS(A) ((A) < 0 ? -(A) : (A))
+#define INDEX2(ROW,COL,NUM_ROWS) ((COL)*(NUM_ROWS)+(ROW))
+
+
+/*======================================================================
+ Main code plus interface to Matlab.
+========================================================================*/
+
+void mexFunction( int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[] )
+{
+ char fname[1000];
+ FILE *fid;
+ char *line;
+ double *feat_val;
+ uint32_t *feat_idx;
+ long nnzf;
+ mxArray *W;
+ long nDims = 0;
+ long nData;
+ long i, j;
+ long nnz = 0;
+ double max_data_norm2 = -mxGetInf();
+ mwIndex *irs, *jcs;
+ double *sr;
+ int verb=0;
+
+ if( nrhs < 1 )
+ mexErrMsgTxt("At least one input argument required.\n"
+ "Synopsis:\n"
+ " [feat,labels] = load_svmlight_format(file_name)\n"
+ " [feat,labels] = load_svmlight_format(file_name,verb)\n"
+ " \n"
+ );
+ if( nrhs >= 2)
+ verb = (int)mxGetScalar(prhs[1]);
+ else
+ verb = 1;
+
+ /* get input arguments */
+ mxGetString(prhs[0], fname, 1000);
+
+ if(verb)
+ mexPrintf("Input file: %s\n", fname);
+
+ fid = fopen(fname, "r");
+ if(fid == NULL) {
+ perror("fopen error: ");
+ mexErrMsgTxt("Cannot open input file.");
+ }
+
+ /**********************************/
+ line = mxCalloc(LIBSLF_MAXLINELEN, sizeof(char));
+ if( line == NULL )
+ mexErrMsgTxt("Not enough memmory to allocate line buffer.");
+
+ feat_idx = mxCalloc(LIBSLF_MAXLINELEN, sizeof(uint32_t));
+ if( feat_idx == NULL )
+ mexErrMsgTxt("Not enough memmory to allocate feat_idx.");
+
+ feat_val = mxCalloc(LIBSLF_MAXLINELEN, sizeof(double));
+ if( feat_val == NULL )
+ mexErrMsgTxt("Not enough memmory to allocate feat_val.");
+
+
+
+ /*********************************************/
+ /* Main code */
+ /*********************************************/
+
+ if(verb)
+ mexPrintf("Analysing input data...");
+
+ double label;
+ int go = 1;
+ long line_cnt = 0;
+
+ while(go) {
+
+ if(fgets(line,LIBSLF_MAXLINELEN, fid) == NULL )
+ {
+ go = 0;
+ if(verb)
+ {
+ if( (line_cnt % 1000) != 0)
+ mexPrintf(" %d", line_cnt);
+ mexPrintf(" EOF.\n");
+ }
+
+ }
+ else
+ {
+ line_cnt ++;
+ nnzf = svmlight_format_parse_line_doubley(line, &label, feat_idx, feat_val);
+
+ if(nnzf == -1)
+ {
+ mexPrintf("Parsing error on line %d .\n", line_cnt);
+ mexErrMsgTxt("Defective input file.");
+ }
+
+ double norm2 = 0;
+ for(j = 0; j < nnzf; j++)
+ norm2 += feat_val[j]*feat_val[j];
+
+ max_data_norm2 = MAX(max_data_norm2,norm2);
+
+ nDims = MAX(nDims,feat_idx[nnzf-1]);
+
+ nnz += nnzf;
+
+ if( (line_cnt % 1000) == 0) {
+ if(verb)
+ {
+ mexPrintf(" %d", line_cnt);
+ fflush(NULL);
+ }
+ }
+ }
+ }
+
+ nData = line_cnt;
+
+ fclose(fid);
+ if(verb)
+ {
+ mexPrintf("Number of examples: %d\n", nData);
+ mexPrintf("Dimensions: %d\n", nDims);
+ mexPrintf("nnz: %d, density: %f%%\n", nnz, 100*(double)nnz/((double)nDims*(double)nData) );
+ mexPrintf("max_i ||x_i||^2: %f\n", max_data_norm2);
+ }
+
+ /*---------------------------------------------*/
+
+
+ mxArray* sp_mat_X = mxCreateSparse(nDims, nData, nnz, mxREAL);
+ if( sp_mat_X == NULL)
+ mexErrMsgTxt("Not enough memory to allocate sp_mat_X");
+ plhs[0] = sp_mat_X;
+
+ plhs[1] = mxCreateDoubleMatrix(nData,1,mxREAL);
+ if( plhs[1] == NULL)
+ mexErrMsgTxt("Not enough memory to allocate vec_y.");
+ double *vec_y = mxGetPr(plhs[1]);
+
+ sr = mxGetPr(sp_mat_X);
+ irs = mxGetIr(sp_mat_X);
+ jcs = mxGetJc(sp_mat_X);
+
+ fid = fopen(fname, "r");
+ if(fid == NULL) {
+ perror("fopen error: ");
+ mexErrMsgTxt("Cannot open input file.");
+ }
+
+ if(verb)
+ mexPrintf("Reading examples...");
+
+ go = 1;
+ line_cnt = 0;
+ long k=0;
+ while(go) {
+ if(fgets(line,LIBSLF_MAXLINELEN, fid) == NULL )
+ {
+ go = 0;
+ if(verb)
+ {
+ if( (line_cnt % 1000) != 0)
+ mexPrintf(" %d", line_cnt);
+ mexPrintf(" EOF.\n");
+ }
+ }
+ else
+ {
+ line_cnt ++;
+ nnzf = svmlight_format_parse_line_doubley(line, &label, feat_idx, feat_val);
+
+ if(nnzf == -1)
+ {
+ mexPrintf("Parsing error on line %d .\n", line_cnt);
+ mexErrMsgTxt("Defective input file.");
+ }
+
+ vec_y[line_cnt-1] = (double)label;
+
+ jcs[line_cnt-1] = k;
+
+ for(j = 0; j < nnzf; j++) {
+ sr[k] = feat_val[j];
+ irs[k] = feat_idx[j]-1;
+ k++;
+ }
+
+ if(verb)
+ {
+ if( (line_cnt % 1000) == 0) {
+ mexPrintf(" %d", line_cnt);
+ fflush(NULL);
+ }
+ }
+ }
+ }
+ jcs[line_cnt] = k;
+
+ plhs[0] = sp_mat_X;
+
+ fclose(fid);
+
+ return;
+}
+
diff --git a/msvmocas b/msvmocas
deleted file mode 100755
index 0d7e283..0000000
Binary files a/msvmocas and /dev/null differ
diff --git a/msvmocas.mexa64 b/msvmocas.mexa64
deleted file mode 100755
index f905495..0000000
Binary files a/msvmocas.mexa64 and /dev/null differ
diff --git a/msvmocas_light.mexa64 b/msvmocas_light.mexa64
deleted file mode 100755
index f7a1d4b..0000000
Binary files a/msvmocas_light.mexa64 and /dev/null differ
diff --git a/ocas_helper.c b/ocas_helper.c
index bf38e65..0e0274e 100644
--- a/ocas_helper.c
+++ b/ocas_helper.c
@@ -78,6 +78,52 @@ static const uint32_t MinimalParallelCutLenght = 100;
/*----------------------------------------------------------------------
+ ----------------------------------------------------------------------*/
+int full_add_nnw_constr(uint32_t idx, uint32_t nSel, void* user_data)
+{
+ full_A[LIBOCAS_INDEX(idx,nSel,nDim)] = 1.0;
+ A0[nSel] = 0.0;
+
+ return( 0 );
+}
+
+
+/*----------------------------------------------------------------------
+ ----------------------------------------------------------------------*/
+int sparse_add_nnw_constr(uint32_t idx, uint32_t nSel, void* user_data)
+{
+ sparse_A.nz_dims[nSel] = 1;
+ sparse_A.index[nSel] = NULL;
+ sparse_A.value[nSel] = NULL;
+ sparse_A.index[nSel] = mxCalloc(1,sizeof(uint32_t));
+ sparse_A.value[nSel] = mxCalloc(1,sizeof(double));
+ if(sparse_A.index[nSel]==NULL || sparse_A.value[nSel]==NULL)
+ {
+ mxFree(sparse_A.index[nSel]);
+ mxFree(sparse_A.value[nSel]);
+ return(-1);
+ }
+
+ sparse_A.index[nSel][0] = idx;
+ sparse_A.value[nSel][0] = 1.0;
+
+ A0[nSel] = 0.0;
+ return( 0 );
+}
+
+
+/*----------------------------------------------------------------------
+ ----------------------------------------------------------------------*/
+void clip_neg_W( uint32_t num_pw_constr, uint32_t *pw_idx, void* user_data )
+{
+ uint32_t i;
+ for(i=0; i< num_pw_constr; i++ ) {
+ W[LIBOCAS_INDEX(pw_idx[i],0,nDim)] = LIBOCAS_MAX(0,W[LIBOCAS_INDEX(pw_idx[i],0,nDim)]);
+ }
+}
+
+
+/*----------------------------------------------------------------------
sq_norm_W = sparse_compute_W( alpha, nSel ) does the following:
oldW = W;
@@ -418,6 +464,7 @@ void full_compute_W( double *sq_norm_W, double *dp_WoldW, double *alpha, uint32_
}
+
/*----------------------------------------------------------------------
parallel_sparse_compute_output( output, user_data ) does the following
diff --git a/ocas_helper.h b/ocas_helper.h
index 8eb14bb..b2cc36a 100644
--- a/ocas_helper.h
+++ b/ocas_helper.h
@@ -57,14 +57,20 @@ extern void ocas_print(ocas_return_value_T value);
extern void ocas_print_null(ocas_return_value_T value);
extern int qsort_data(double* value, double* data, uint32_t size);
-/* helper functions for two-class SVM solver */
+/** helper functions for two-class SVM solver **/
extern double update_W( double t, void* user_data );
/* dense cutting plane buffer*/
extern void full_compute_W( double *sq_norm_W, double *dp_WoldW, double *alpha, uint32_t nSel, void* user_data );
/* sparse cutting plane bugger */
extern void sparse_compute_W( double *sq_norm_W, double *dp_WoldW, double *alpha, uint32_t nSel, void* user_data );
-/* helper functions for multi-class SVM solver */
+/** helper function for two-class SVM solver with additional positivity contraints on W **/
+int full_add_nnw_constr(uint32_t idx, uint32_t nSel, void* user_data);
+void clip_neg_W( uint32_t num_pw_constr, uint32_t *pw_idx, void* user_data );
+int sparse_add_nnw_constr(uint32_t idx, uint32_t nSel, void* user_data);
+
+
+/** helper functions for multi-class SVM solver **/
extern double msvm_update_W( double t, void* user_data );
/* dense cutting plane buffer*/
extern void msvm_full_compute_W( double *sq_norm_W, double *dp_WoldW, double *alpha, uint32_t nSel, void* user_data );
diff --git a/svmocas b/svmocas
deleted file mode 100755
index df8c752..0000000
Binary files a/svmocas and /dev/null differ
diff --git a/svmocas.mexa64 b/svmocas.mexa64
deleted file mode 100755
index d854b06..0000000
Binary files a/svmocas.mexa64 and /dev/null differ
diff --git a/svmocas_bool_mex.c b/svmocas_bool_mex.c
new file mode 100644
index 0000000..3231781
--- /dev/null
+++ b/svmocas_bool_mex.c
@@ -0,0 +1,314 @@
+/*=================================================================================
+ * svmocas_bool_mex.c: Matlab MEX interface for OCAS solver training the
+ * linear SVM classifiers from boolean features.
+ *
+ * Synopsis:
+ * [W,W0,stat] = svmocas_bool(X,X0,y,C,Method,TolRel,TolAbs,QPBound,BufSize,nData,MaxTime,verb)
+ *
+ * Copyright (C) 2011 Vojtech Franc, xfrancv at cmp.felk.cvut.cz
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public
+ * License as published by the Free Software Foundation;
+ *======================================================================================*/
+
+#include <stdio.h>
+#include <string.h>
+#include <stdint.h>
+#include <mex.h>
+
+#include "libocas.h"
+#include "ocas_helper.h"
+#include "features_bool.h"
+
+#define DEFAULT_METHOD 1
+#define DEFAULT_TOLREL 0.01
+#define DEFAULT_TOLABS 0.0
+#define DEFAULT_QPVALUE 0.0
+#define DEFAULT_BUFSIZE 2000
+#define DEFAULT_MAXTIME mxGetInf()
+#define DEFAULT_VERB 1
+
+/*======================================================================
+ Main code plus interface to Matlab.
+========================================================================*/
+
+void mexFunction( int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[] )
+{
+ double C, TolRel, TolAbs, QPBound, trn_err, MaxTime;
+ double *vec_C;
+ int8_t *ptr;
+ uint32_t num_of_Cs;
+ uint32_t i, j, BufSize;
+ uint16_t Method;
+ int verb;
+ ocas_return_value_T ocas;
+
+ /* timing variables */
+ double init_time;
+ double total_time;
+
+ total_time = get_time();
+ init_time = total_time;
+
+ if(nrhs < 4 || nrhs > 12)
+ mexErrMsgTxt("Improper number of input arguments.\n\n"
+ "SVMOCAS_BOOL solver for training two-class linear SVM classifiers\n\n"
+ "Synopsis:\n"
+ " [W,W0,stat] = svmocas_bool(X,X0,y,C,Method,TolRel,TolAbs,QPBound,"
+ "BufSize,nExamples,MaxTime) \n\n"
+ "Input: \n"
+ " X [nDim x nExamples] \n"
+ " X0 [1 x 1 (double)] constant feature added to all examples\n"
+ " y [nExamples x 1 (double)] labels of the examples (+1/-1)\n"
+ " C [1x1] or [nExamples x 1] regularization constant(s) \n"
+ " Method [1x1 (double)] 0 for BMRM; 1 for OCAS \n"
+ " TolRel [1x1 (double)]\n"
+ " TolAbs [1x1 (double)]\n"
+ " QPBound [1x1 (double)]\n"
+ " BufSize [1x1 (double)]\n"
+ " nExamples [1x1 (double) number of examples to use; "
+ "(inf means use all examples)\n"
+ " MaxTime [1x1 (double)]\n"
+ " verb [1x1 (bouble)]\n\n"
+ "Output:\n"
+ " W [nDim x 1] Parameter vector\n"
+ " W0 [1x1] Bias term\n"
+ " stat [struct] \n");
+
+ data_X = (mxArray*)prhs[0];
+
+ /* check the input matrix X */
+ if( !mxIsUint8(data_X) || mxIsSparse(data_X) )
+ mexErrMsgTxt("The features X must be a dense matrix of uint8.");
+
+ if (mxGetNumberOfDimensions(data_X) != 2)
+ mexErrMsgTxt("Input argument X must be two dimensional");
+
+ if(nrhs >= 12)
+ verb = (int)mxGetScalar(prhs[11]);
+ else
+ verb = DEFAULT_VERB;
+
+ X0 = (double)mxGetScalar(prhs[1]);
+
+ if( !mxIsDouble(prhs[2]))
+ mexErrMsgTxt("The labels y must be a dense vector of doubles.");
+
+ data_y = (double*)mxGetPr(prhs[2]);
+
+ if(LIBOCAS_MAX(mxGetM(prhs[2]),mxGetN(prhs[2])) != mxGetN(prhs[0]))
+ mexErrMsgTxt("Length of vector y must equl to the number of columns of matrix X.");
+
+ /********************************************/
+ /* zmenit */
+ /********************************************/
+ nDim = mxGetM(prhs[0])*8;
+/* nDim = mxGetM(prhs[0]);*/
+
+ if(verb)
+ {
+ mexPrintf("Input data statistics:\n"
+ " # of examples : %d\n"
+ " dimensionality : %d\n",
+ mxGetN(data_X), nDim);
+ }
+
+
+ num_of_Cs = LIBOCAS_MAX(mxGetN(prhs[3]),mxGetM(prhs[3]));
+ if(num_of_Cs == 1)
+ {
+ C = (double)mxGetScalar(prhs[3]);
+ }
+ else
+ {
+ vec_C = (double*)mxGetPr(prhs[3]);
+ }
+
+ if(nrhs >= 5)
+ Method = (uint32_t)mxGetScalar(prhs[4]);
+ else
+ Method = DEFAULT_METHOD;
+
+ if(nrhs >= 6)
+ TolRel = (double)mxGetScalar(prhs[5]);
+ else
+ TolRel = DEFAULT_TOLREL;
+
+ if(nrhs >= 7)
+ TolAbs = (double)mxGetScalar(prhs[6]);
+ else
+ TolAbs = DEFAULT_TOLABS;
+
+ if(nrhs >= 8)
+ QPBound = (double)mxGetScalar(prhs[7]);
+ else
+ QPBound = DEFAULT_QPVALUE;
+
+ if(nrhs >= 9)
+ BufSize = (uint32_t)mxGetScalar(prhs[8]);
+ else
+ BufSize = DEFAULT_BUFSIZE;
+
+ if(nrhs >= 10 && mxIsInf(mxGetScalar(prhs[9])) == false)
+ nData = (uint32_t)mxGetScalar(prhs[9]);
+ else
+ nData = mxGetN(data_X);
+
+ if(nData < 1 || nData > mxGetN(prhs[0]))
+ mexErrMsgTxt("Improper value of argument nData.");
+
+ if(num_of_Cs > 1 && num_of_Cs < nData)
+ mexErrMsgTxt("Length of the vector C less than the number of examples.");
+
+ if(nrhs >= 11)
+ MaxTime = (double)mxGetScalar(prhs[10]);
+ else
+ MaxTime = DEFAULT_MAXTIME;
+
+
+
+ /*----------------------------------------------------------------
+ Print setting
+ -------------------------------------------------------------------*/
+ if(verb)
+ {
+ mexPrintf("Setting:\n");
+
+ if( num_of_Cs == 1)
+ mexPrintf(" C : %f\n", C);
+ else
+ mexPrintf(" C : different for each example\n");
+
+ mexPrintf(" bias : %.0f\n"
+ " # of examples : %d\n"
+ " solver : %d\n"
+ " cache size : %d\n"
+ " TolAbs : %f\n"
+ " TolRel : %f\n"
+ " QPValue : %f\n"
+ " MaxTime : %f [s]\n"
+ " verb : %d\n",
+ X0, nData, Method,BufSize,TolAbs,TolRel, QPBound, MaxTime, verb);
+ }
+
+ /* learned weight vector */
+ plhs[0] = (mxArray*)mxCreateDoubleMatrix(nDim,1,mxREAL);
+ W = (double*)mxGetPr(plhs[0]);
+ if(W == NULL) mexErrMsgTxt("Not enough memory for vector W.");
+
+ oldW = (double*)mxCalloc(nDim,sizeof(double));
+ if(oldW == NULL) mexErrMsgTxt("Not enough memory for vector oldW.");
+
+ W0 = 0;
+ oldW0 = 0;
+
+ A0 = mxCalloc(BufSize,sizeof(A0[0]));
+ if(A0 == NULL) mexErrMsgTxt("Not enough memory for vector A0.");
+
+ /* allocate buffer for computing cutting plane */
+ new_a = (double*)mxCalloc(nDim,sizeof(double));
+ if(new_a == NULL)
+ mexErrMsgTxt("Not enough memory for auxiliary cutting plane buffer new_a.");
+
+ if(num_of_Cs > 1)
+ {
+ for(i=0; i < nData; i++)
+ data_y[i] = data_y[i]*vec_C[i];
+ }
+
+ /* init cutting plane buffer */
+ full_A = mxCalloc(BufSize*nDim,sizeof(double));
+ if( full_A == NULL )
+ mexErrMsgTxt("Not enough memory for cutting plane buffer full_A.");
+
+
+ /* select function to print progress info */
+ void (*print_function)(ocas_return_value_T);
+ if(verb)
+ {
+ mexPrintf("Starting optimization:\n");
+ print_function = &ocas_print;
+ }
+ else
+ {
+ print_function = &ocas_print_null;
+ }
+
+ init_time=get_time()-init_time;
+
+ if(num_of_Cs == 1)
+ ocas = svm_ocas_solver( C, nData, TolRel, TolAbs, QPBound, MaxTime,BufSize, Method,
+ &full_compute_W, &update_W, &full_bool_add_new_cut,
+ &full_bool_compute_output, &qsort_data, print_function, 0);
+ else
+ ocas = svm_ocas_solver_difC( vec_C, nData, TolRel, TolAbs, QPBound, MaxTime,BufSize, Method,
+ &full_compute_W, &update_W, &full_bool_add_new_cut,
+ &full_bool_compute_output, &qsort_data, print_function, 0);
+
+ total_time=get_time()-total_time;
+
+ if(verb)
+ {
+ mexPrintf("Stopping condition: ");
+ switch( ocas.exitflag )
+ {
+ case 1: mexPrintf("1-Q_D/Q_P <= TolRel(=%f) satisfied.\n", TolRel); break;
+ case 2: mexPrintf("Q_P-Q_D <= TolAbs(=%f) satisfied.\n", TolAbs); break;
+ case 3: mexPrintf("Q_P <= QPBound(=%f) satisfied.\n", QPBound); break;
+ case 4: mexPrintf("Optimization time (=%f) >= MaxTime(=%f).\n", ocas.ocas_time, MaxTime);
+ break;
+ case -1: mexPrintf("Has not converged!\n" ); break;
+ case -2: mexPrintf("Not enough memory for the solver.\n" ); break;
+ }
+
+ mexPrintf("Timing statistics:\n"
+ " init_time : %f[s]\n"
+ " qp_solver_time : %f[s]\n"
+ " sort_time : %f[s]\n"
+ " output_time : %f[s]\n"
+ " add_time : %f[s]\n"
+ " w_time : %f[s]\n"
+ " print_time : %f[s]\n"
+ " ocas_time : %f[s]\n"
+ " total_time : %f[s]\n",
+ init_time, ocas.qp_solver_time, ocas.sort_time, ocas.output_time,
+ ocas.add_time, ocas.w_time, ocas.print_time, ocas.ocas_time, total_time);
+
+ mexPrintf("Training error: %.4f%%\n", 100*(double)ocas.trn_err/(double)nData);
+ }
+
+ if(num_of_Cs > 1)
+ {
+ for(i=0; i < nData; i++)
+ data_y[i] = data_y[i]/vec_C[i];
+ }
+
+ plhs[1] = mxCreateDoubleScalar( W0 );
+
+ const char *field_names[] = {"nTrnErrors","Q_P","Q_D","nIter","nCutPlanes","exitflag",
+ "init_time","output_time","sort_time","qp_solver_time","add_time",
+ "w_time","print_time","ocas_time","total_time"};
+ mwSize dims[2] = {1,1};
+
+ plhs[2] = mxCreateStructArray(2, dims, (sizeof(field_names)/sizeof(*field_names)), field_names);
+
+ mxSetField(plhs[2],0,"nIter",mxCreateDoubleScalar((double)ocas.nIter));
+ mxSetField(plhs[2],0,"nCutPlanes",mxCreateDoubleScalar((double)ocas.nCutPlanes));
+ mxSetField(plhs[2],0,"nTrnErrors",mxCreateDoubleScalar(ocas.trn_err));
+ mxSetField(plhs[2],0,"Q_P",mxCreateDoubleScalar(ocas.Q_P));
+ mxSetField(plhs[2],0,"Q_D",mxCreateDoubleScalar(ocas.Q_D));
+ mxSetField(plhs[2],0,"init_time",mxCreateDoubleScalar(init_time));
+ mxSetField(plhs[2],0,"output_time",mxCreateDoubleScalar(ocas.output_time));
+ mxSetField(plhs[2],0,"sort_time",mxCreateDoubleScalar(ocas.sort_time));
+ mxSetField(plhs[2],0,"qp_solver_time",mxCreateDoubleScalar(ocas.qp_solver_time));
+ mxSetField(plhs[2],0,"add_time",mxCreateDoubleScalar(ocas.add_time));
+ mxSetField(plhs[2],0,"w_time",mxCreateDoubleScalar(ocas.w_time));
+ mxSetField(plhs[2],0,"print_time",mxCreateDoubleScalar(ocas.print_time));
+ mxSetField(plhs[2],0,"ocas_time",mxCreateDoubleScalar(ocas.ocas_time));
+ mxSetField(plhs[2],0,"total_time",mxCreateDoubleScalar(total_time));
+ mxSetField(plhs[2],0,"exitflag",mxCreateDoubleScalar((double)ocas.exitflag));
+
+ return;
+}
+
diff --git a/svmocas_bool_test.m b/svmocas_bool_test.m
new file mode 100644
index 0000000..7ef1f08
--- /dev/null
+++ b/svmocas_bool_test.m
@@ -0,0 +1,16 @@
+%
+clc;
+C=100;
+
+X=int8([0 1 0 1;0 0 1 1;zeros(6,4);0 1 0 1;0 0 1 1;zeros(6,4)]);
+y = [+1 +1 -1 -1];
+[W,W0,stat] = svmocas(X,1,y,C);
+
+%%%
+
+Xbool = uint8([0 1 2 3; 0 1 2 3]);
+y = [+1 +1 -1 -1];
+[Wbool,W0bool,stat] = svmocas_bool(Xbool,1,y,C);
+
+W-Wbool
+W0-W0bool
diff --git a/svmocas_lbp.mexa64 b/svmocas_lbp.mexa64
deleted file mode 100755
index e39e4ae..0000000
Binary files a/svmocas_lbp.mexa64 and /dev/null differ
diff --git a/svmocas_light.mexa64 b/svmocas_light.mexa64
deleted file mode 100755
index f62cee8..0000000
Binary files a/svmocas_light.mexa64 and /dev/null differ
diff --git a/svmocas_nnw.m b/svmocas_nnw.m
new file mode 100644
index 0000000..ad7bec3
--- /dev/null
+++ b/svmocas_nnw.m
@@ -0,0 +1,74 @@
+% SVMOCAS_NNW training linear two-class SVM classifiers with non-negative weights
+%
+% Synopsis:
+% [W,W0,stat] = svmocas_nnw(X,X0,y,C,nnw,Method,TolRel,TolAbs,QPBound,BufSize,nExamples,MaxTime)
+%
+% Desription:
+% This function trains linear SVM classifier with non-negatie weight by solving
+%
+% W,W0 = argmin 0.5*(w'*w+w0^2) + C*sum max( 0, 1-y(i)*(w'*X(:,i)+w0*X0) )
+% w,w0 i=1:nExamples
+%
+% subject to w(nnw) >= 0
+%
+% The standard SVM is obtained if nnw is empty, however, in this case use
+% the function SVMOCAS instead.
+%
+% The training examples are passed in the matrix X where each column is a single
+% feature vector. The matrix X can be one of the following types
+% dense double (default in Matlab)
+% sparse double
+% dense single
+% dense int8
+% The labels of the training examples must be given as dense vector y whose length
+% equals to the number of columns of the matrix X.
+%
+% Reference:
+% V. Franc, S. Sonnenburg. Optimized Cutting Plane Algorithm for Large-scale Risk
+% Minimization. Journal of Machine Learning Research. 2009
+% http://jmlr.csail.mit.edu/papers/volume10/franc09a/franc09a.pdf
+%
+% Input:
+% X [nDim x nExamples] training inputs (dense double, sparse double, dense single,
+% dense int8 matrix).
+% X0 [1x1] constant coordinate (implicitly) added to all examples;
+% this allows training biased decision rule.
+% y [nExamples x 1] labels (+1/-1).
+% C [1x1] or [nExamples x 1] C [1x1] is a positive (>0) regularization constant.
+% nnz [P x 1] Indices of weights of W which must be non-negative.
+% Method [1x1] 0..cutting plane; 1..OCAS (default 1)
+% TolRel [1x1] halts if Q_P-Q_D <= abs(Q_P)*TolRel (default 0.01)
+% TolAbs [1x1] halts if Q_P-Q_D <= TolAbs (default 0)
+% QPValue [1x1] halts if Q_P <= QPBpound (default 0)
+% BufSize [1x1] Initial size of active constrains buffer (default 2000)
+% nExaples [1x1] Number of examples used for training; must be >0 and <= size(X,2).
+% If nExamples = inf then nExamples is set to size(X,2).
+% MaxTime [1x1] halts if time used by solver (data loading time is not counted) exceeds
+% MaxTime given in seconds. Use MaxTime=inf (default) to switch off this stopping condition.
+%
+% Output:
+% W [nDim x 1] Paramater vectors of decision rule sign(W'*X+W0)
+% W0 [1x1] Bias term of the decision rule.
+% stat [struct] Optimizer statistics (field names are self-explaining).
+%
+% Example:
+% % train standard linear SVM classifier and then SVM classifier whose
+% % first weight is non-negative.
+%
+% libocasPath = fileparts(which('svmocas'));
+% trn = load([libocasPath '/data/riply_trn.mat'],'X','y');
+% svmC = 1;
+% [W,W0,stat] = svmocas(trn.X,1,trn.y,svmC);
+% [W_nn,W0_nn,stat] = svmocas_nnw(trn.X,1,trn.y,svmC,[1]);
+%
+% W
+% W_nn
+%
+
+%
+% Copyright (C) 2008, 2009, 2012 Vojtech Franc, xfrancv at cmp.felk.cvut.cz
+% Soeren Sonnenburg, soeren.sonnenburg at first.fraunhofer.de
+%
+% This program is free software; you can redistribute it and/or
+% modify it under the terms of the GNU General Public
+% License as published by the Free Software Foundation;
diff --git a/svmocas_nnw_mex.c b/svmocas_nnw_mex.c
new file mode 100644
index 0000000..bb65020
--- /dev/null
+++ b/svmocas_nnw_mex.c
@@ -0,0 +1,456 @@
+/*======================================================================================
+ * svmocas_mex.c: Matlab MEX interface to OCAS solver for training two-class
+ * linear SVM classifier
+ *
+ * Synopsis:
+ * [W,W0,stat] = svmocas_nnw(X,X0,y,C,nnw,Method,TolRel,TolAbs,QPBound,BufSize,
+ * nData,MaxTime,verb)
+ *
+ * See svmocas.m for more help.
+ *
+ * Copyright (C) 2008, 2009 Vojtech Franc, xfrancv at cmp.felk.cvut.cz
+ * Soeren Sonnenburg, soeren.sonnenburg at first.fraunhofer.de
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public
+ * License as published by the Free Software Foundation;
+ *=====================================================================================*/
+
+#include <stdio.h>
+#include <string.h>
+#include <stdint.h>
+#include <mex.h>
+
+#include "libocas.h"
+#include "ocas_helper.h"
+#include "features_int8.h"
+#include "features_double.h"
+#include "features_single.h"
+
+#define DEFAULT_METHOD 1
+#define DEFAULT_TOLREL 0.01
+#define DEFAULT_TOLABS 0.0
+#define DEFAULT_QPVALUE 0.0
+#define DEFAULT_BUFSIZE 2000
+#define DEFAULT_MAXTIME mxGetInf()
+#define DEFAULT_VERB 1
+
+/*======================================================================
+ Main code plus interface to Matlab.
+========================================================================*/
+
+void mexFunction( int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[] )
+{
+ double C, TolRel, TolAbs, QPBound, trn_err, MaxTime;
+ uint32_t i, j, BufSize;
+ uint16_t Method;
+ int verb;
+ double *ptr_double;
+ float *ptr_single;
+ int8_t *ptr_int8;
+ int32_t num_nnw, *nnw_idx=NULL;
+ ocas_return_value_T ocas;
+
+ /* timing variables */
+ double init_time;
+ double total_time;
+
+ total_time = get_time();
+ init_time = total_time;
+
+ if(nrhs < 5 || nrhs > 13)
+ mexErrMsgTxt("Improper number of input arguments.\n\n"
+ "SVMOCAS_NNW training two-class linear SVM classifiers with non-negative weights\n\n"
+ "Synopsis:\n"
+ " [W,W0,stat] = svmocas_nnw(X,X0,y,C,nnw,Method,TolRel,TolAbs,QPBound,"
+ "BufSize,nExamples,MaxTime,verb) \n\n"
+ "Input: \n"
+ " X [nDim x nExamples] training inputs (dense double or sparse double or dense single or dense int8)\n"
+ " X0 [1 x 1 (double)] constant feature added to all examples\n"
+ " y [nExamples x 1 (double)] labels of the examples (+1/-1)\n"
+ " C [1x1] a positive regularization constant \n"
+ " nnw [p x 1 (double)] a vector of indices of weights W which should be non-negative (p <= nDim) \n"
+ " Method [1x1 (double)] 0 for BMRM; 1 for OCAS \n"
+ " TolRel [1x1 (double)] stop if F(W)-F(W_optimal) <= abs(F(W))*TolRel\n"
+ " TolAbs [1x1 (double)] stop if F(W)-F(W_optimal) <= TolAbs\n"
+ " QPBound [1x1 (double)] stop is F(W) <= QPBound\n"
+ " BufSize [1x1 (double)] size of cutting plane buffer\n"
+ " nExamples [1x1 (double) number of examples to use; (inf means use all examples)\n"
+ " MaxTime [1x1 (double)] stop if processing time exceeds MaxTime given in seconds\n"
+ " verb [1x1 (bouble)] if 1 then it prints progress and other info\n\n"
+ "Output:\n"
+ " W [nDim x 1] Parameter vector\n"
+ " W0 [1x1] Bias term\n"
+ " stat [struct] Optimization statistics (time, num of iterations etc)\n");
+
+
+ if(nrhs >= 13)
+ verb = (int)mxGetScalar(prhs[11]);
+ else
+ verb = DEFAULT_VERB;
+
+ /* 1st input argument: training feature vectors */
+ data_X = (mxArray*)prhs[0];
+ if( (mxGetNumberOfDimensions(data_X) != 2) ||
+ !( ( mxIsDouble(data_X) && mxIsSparse(data_X) ) ||
+ ( mxIsDouble(data_X) && !mxIsSparse(data_X) ) ||
+ ( mxIsSingle(data_X) && !mxIsSparse(data_X) ) ||
+ ( mxIsInt8(data_X) && !mxIsSparse(data_X) ) ))
+ {
+ mexErrMsgTxt("The first input argument must be two dimensional matrix of the following type:\n"
+ "dense double or sparse double or dense single or dense int8 matrix.\n");
+ }
+
+ /* 2nd input argument: constant coordinate added to feature vectors */
+ X0 = (double)mxGetScalar(prhs[1]);
+
+ /*3rd input argument: vector of labels */
+ if( !mxIsDouble(prhs[2]) || mxIsSparse(prhs[2]) )
+ mexErrMsgTxt("The third input argument must be dense vector of doubles.");
+ data_y = (double*)mxGetPr(prhs[2]);
+
+ if(LIBOCAS_MAX(mxGetM(prhs[2]),mxGetN(prhs[2])) != mxGetN(prhs[0]))
+ mexErrMsgTxt("Length of vector y (3rd input argument) must be equal to the number of columns of matrix X (1st input argument).");
+
+ nDim = mxGetM(prhs[0]);
+
+ if(verb)
+ {
+ mexPrintf("Input data statistics:\n"
+ " # of examples : %d\n"
+ " dimensionality : %d\n",
+ mxGetN(data_X), nDim);
+
+ if( mxIsSparse(data_X)== true )
+ mexPrintf(" sparse features (density=%.2f%%) ",
+ 100.0*(double)mxGetNzmax(data_X)/((double)nDim*(double)(mxGetN(data_X))));
+ else
+ mexPrintf(" dense features ");
+ if( mxIsDouble(data_X) )
+ mexPrintf("in double precision\n");
+ if( mxIsSingle(data_X) )
+ mexPrintf("in single precision\n");
+ if( mxIsInt8(data_X) )
+ mexPrintf("represented as int8\n");
+ }
+
+ if(LIBOCAS_MAX(mxGetN(prhs[3]),mxGetM(prhs[3])) == 1)
+ C = (double)mxGetScalar(prhs[3]);
+ else
+ mexErrMsgTxt("The fourth argumetn C must be a positive scalar.");
+
+ num_nnw = LIBOCAS_MAX(mxGetN(prhs[4]),mxGetM(prhs[4]));
+ if( num_nnw < 1 )
+ mexErrMsgTxt("The fifth argument nnw must be a non-empty vector.");
+
+ ptr_double = mxGetPr(prhs[4]);
+ nnw_idx = (uint32_t*)mxCalloc(num_nnw,sizeof(uint32_t));
+ if(nnw_idx == NULL) mexErrMsgTxt("Not enough memory for vector nnw_idx.");
+ for( i = 0; i < num_nnw; i++) {
+ if( ptr_double[i] < 1 || ptr_double[i] > nDim)
+ mexErrMsgTxt("The fifth argument nnw containts improper values.");
+
+ nnw_idx[i] = (uint32_t)ptr_double[i]-1;
+ }
+
+
+
+ if(nrhs >= 6)
+ Method = (uint32_t)mxGetScalar(prhs[5]);
+ else
+ Method = DEFAULT_METHOD;
+
+ if(nrhs >= 7)
+ TolRel = (double)mxGetScalar(prhs[6]);
+ else
+ TolRel = DEFAULT_TOLREL;
+
+ if(nrhs >= 8)
+ TolAbs = (double)mxGetScalar(prhs[7]);
+ else
+ TolAbs = DEFAULT_TOLABS;
+
+ if(nrhs >= 9)
+ QPBound = (double)mxGetScalar(prhs[8]);
+ else
+ QPBound = DEFAULT_QPVALUE;
+
+ if(nrhs >= 10)
+ BufSize = (uint32_t)mxGetScalar(prhs[9]);
+ else
+ BufSize = DEFAULT_BUFSIZE;
+
+ if(nrhs >= 11 && mxIsInf(mxGetScalar(prhs[10])) == false)
+ nData = (uint32_t)mxGetScalar(prhs[10]);
+ else
+ nData = mxGetN(data_X);
+
+ if(nData < 1 || nData > mxGetN(prhs[0]))
+ mexErrMsgTxt("Improper value of argument nData.");
+
+ if(nrhs >= 12)
+ MaxTime = (double)mxGetScalar(prhs[11]);
+ else
+ MaxTime = DEFAULT_MAXTIME;
+
+
+ /*----------------------------------------------------------------
+ Print setting
+ -------------------------------------------------------------------*/
+ if(verb)
+ {
+ mexPrintf("Setting:\n");
+
+ mexPrintf(" C : %f\n", C);
+ mexPrintf(" bias : %.0f\n"
+ " # of examples : %d\n"
+ " solver : %d\n"
+ " cache size : %d\n"
+ " TolAbs : %f\n"
+ " TolRel : %f\n"
+ " QPValue : %f\n"
+ " MaxTime : %f [s]\n"
+ " verb : %d\n",
+ X0, nData, Method,BufSize,TolAbs,TolRel, QPBound, MaxTime, verb);
+ }
+
+ /* learned weight vector */
+ plhs[0] = (mxArray*)mxCreateDoubleMatrix(nDim,1,mxREAL);
+ W = (double*)mxGetPr(plhs[0]);
+ if(W == NULL) mexErrMsgTxt("Not enough memory for vector W.");
+
+ oldW = (double*)mxCalloc(nDim,sizeof(double));
+ if(oldW == NULL) mexErrMsgTxt("Not enough memory for vector oldW.");
+
+ W0 = 0;
+ oldW0 = 0;
+
+ A0 = mxCalloc(BufSize,sizeof(A0[0]));
+ if(A0 == NULL) mexErrMsgTxt("Not enough memory for vector A0.");
+
+ /* allocate buffer for computing cutting plane */
+ new_a = (double*)mxCalloc(nDim,sizeof(double));
+ if(new_a == NULL)
+ mexErrMsgTxt("Not enough memory for auxciliary cutting plane buffer new_a.");
+
+ /* select function to print progress info */
+ void (*print_function)(ocas_return_value_T);
+ if(verb)
+ {
+ mexPrintf("Starting optimization:\n");
+ print_function = &ocas_print;
+ }
+ else
+ {
+ print_function = &ocas_print_null;
+ }
+
+ if( mxIsSparse(data_X)== true )
+ {
+
+ /* for i=1:nData, X(:,i) = X(:,i)*y(i); end*/
+ for(i=0; i < nData; i++)
+ mul_sparse_col(data_y[i], data_X, i);
+
+
+ /* init cutting plane buffer */
+ sparse_A.nz_dims = mxCalloc(BufSize,sizeof(uint32_t));
+ sparse_A.index = mxCalloc(BufSize,sizeof(sparse_A.index[0]));
+ sparse_A.value = mxCalloc(BufSize,sizeof(sparse_A.value[0]));
+ if(sparse_A.nz_dims == NULL || sparse_A.index == NULL || sparse_A.value == NULL)
+ mexErrMsgTxt("Not enough memory for cutting plane buffer sparse_A.");
+
+ init_time=get_time()-init_time;
+
+
+ /* ocas = svm_ocas_solver( C, nData, TolRel, TolAbs, QPBound, MaxTime,BufSize, Method,
+ &sparse_compute_W, &update_W, &sparse_add_new_cut,
+ &sparse_compute_output, &qsort_data,
+ print_function, 0);*/
+
+ ocas = svm_ocas_solver_nnw( C, nData, num_nnw, nnw_idx, TolRel, TolAbs, QPBound, MaxTime,
+ BufSize, Method,
+ &sparse_add_nnw_constr, &clip_neg_W, &sparse_compute_W,
+ &update_W, &sparse_add_new_cut,
+ &sparse_compute_output, &qsort_data, print_function, 0);
+
+
+
+ }
+ else
+ {
+
+ int (*add_new_cut)(double*, uint32_t*, uint32_t, uint32_t, void*);
+ int (*compute_output)( double*, void* );
+
+ /* features in double precision */
+ if( mxIsDouble(data_X) )
+ {
+ ptr_double = mxGetPr(data_X);
+ for(i=0; i < nData; i++) {
+ for(j=0; j < nDim; j++ ) {
+ ptr_double[LIBOCAS_INDEX(j,i,nDim)] = ptr_double[LIBOCAS_INDEX(j,i,nDim)]*data_y[i];
+ }
+ }
+
+ add_new_cut = &full_add_new_cut;
+ compute_output = &full_compute_output;
+ }
+
+ /* features in single precision */
+ if( mxIsSingle(data_X) )
+ {
+ mexPrintf("Multiplying data with labels for single precision.n");
+ ptr_single = (float*)mxGetPr(data_X);
+ for(i=0; i < nData; i++) {
+ for(j=0; j < nDim; j++ ) {
+ ptr_single[LIBOCAS_INDEX(j,i,nDim)] = ptr_single[LIBOCAS_INDEX(j,i,nDim)]*(float)data_y[i];
+ }
+ }
+
+ add_new_cut = &full_single_add_new_cut;
+ compute_output = &full_single_compute_output;
+ }
+
+ /* features in int8 */
+ if( mxIsInt8(data_X) )
+ {
+ ptr_int8 = (int8_t*)mxGetPr(data_X);
+ for(i=0; i < nData; i++) {
+ for(j=0; j < nDim; j++ ) {
+ ptr_int8[LIBOCAS_INDEX(j,i,nDim)] = ptr_int8[LIBOCAS_INDEX(j,i,nDim)]/(int8_t)data_y[i];
+ }
+ }
+
+ add_new_cut = &full_int8_add_new_cut;
+ compute_output = &full_int8_compute_output;
+ }
+
+ /* init cutting plane buffer */
+ full_A = mxCalloc(BufSize*nDim,sizeof(double));
+ if( full_A == NULL )
+ mexErrMsgTxt("Not enough memory for cutting plane buffer full_A.");
+
+ init_time=get_time()-init_time;
+
+ /* ocas = svm_ocas_solver( C, nData, TolRel, TolAbs, QPBound, MaxTime,BufSize, Method,
+ &full_compute_W, &update_W, add_new_cut,
+ compute_output, &qsort_data, print_function, 0);
+ */
+
+ ocas = svm_ocas_solver_nnw( C, nData, num_nnw, nnw_idx, TolRel, TolAbs, QPBound, MaxTime,
+ BufSize, Method,
+ &full_add_nnw_constr, &clip_neg_W, &full_compute_W,
+ &update_W, add_new_cut,
+ compute_output, &qsort_data, print_function, 0);
+ }
+
+ total_time=get_time()-total_time;
+
+ if(verb)
+ {
+ mexPrintf("Stopping condition: ");
+ switch( ocas.exitflag )
+ {
+ case 1: mexPrintf("1-Q_D/Q_P <= TolRel(=%f) satisfied.\n", TolRel); break;
+ case 2: mexPrintf("Q_P-Q_D <= TolAbs(=%f) satisfied.\n", TolAbs); break;
+ case 3: mexPrintf("Q_P <= QPBound(=%f) satisfied.\n", QPBound); break;
+ case 4: mexPrintf("Optimization time (=%f) >= MaxTime(=%f).\n",
+ ocas.ocas_time, MaxTime); break;
+ case -1: mexPrintf("Has not converged!\n" ); break;
+ case -2: mexPrintf("Not enough memory for the solver.\n" ); break;
+ }
+
+ mexPrintf("Timing statistics:\n"
+ " init_time : %f[s]\n"
+ " qp_solver_time : %f[s]\n"
+ " sort_time : %f[s]\n"
+ " output_time : %f[s]\n"
+ " add_time : %f[s]\n"
+ " w_time : %f[s]\n"
+ " print_time : %f[s]\n"
+ " ocas_time : %f[s]\n"
+ " total_time : %f[s]\n",
+ init_time, ocas.qp_solver_time, ocas.sort_time, ocas.output_time,
+ ocas.add_time, ocas.w_time, ocas.print_time, ocas.ocas_time, total_time);
+
+ mexPrintf("Training error: %.4f%%\n", 100*(double)ocas.trn_err/(double)nData);
+ }
+
+ /* multiply data by labels as it was at the begining */
+ if( mxIsSparse(data_X)== true )
+ {
+ /* for i=1:nData, X(:,i) = X(:,i)*y(i); end*/
+ for(i=0; i < nData; i++)
+ {
+ mul_sparse_col(1/data_y[i], data_X, i);
+ }
+ }
+ else
+ {
+
+ /* features in double precision */
+ if( mxIsDouble(data_X) )
+ {
+ ptr_double = mxGetPr(data_X);
+ for(i=0; i < nData; i++) {
+ for(j=0; j < nDim; j++ ) {
+ ptr_double[LIBOCAS_INDEX(j,i,nDim)] = ptr_double[LIBOCAS_INDEX(j,i,nDim)]/data_y[i];
+ }
+ }
+ }
+
+ /* features in single precision */
+ if( mxIsSingle(data_X) )
+ {
+ ptr_single = (float*)mxGetPr(data_X);
+ for(i=0; i < nData; i++) {
+ for(j=0; j < nDim; j++ ) {
+ ptr_single[LIBOCAS_INDEX(j,i,nDim)] = ptr_single[LIBOCAS_INDEX(j,i,nDim)]/data_y[i];
+ }
+ }
+ }
+
+ /* features in int8 */
+ if( mxIsInt8(data_X) )
+ {
+ ptr_int8 = (int8_t*)mxGetPr(data_X);
+ for(i=0; i < nData; i++) {
+ for(j=0; j < nDim; j++ ) {
+ ptr_int8[LIBOCAS_INDEX(j,i,nDim)] = ptr_int8[LIBOCAS_INDEX(j,i,nDim)]/data_y[i];
+ }
+ }
+ }
+ }
+
+ /* create output variables */
+ plhs[1] = mxCreateDoubleScalar( W0 );
+
+ const char *field_names[] = {"nTrnErrors","Q_P","Q_D","nIter","nCutPlanes","exitflag",
+ "init_time","output_time","sort_time",
+ "qp_solver_time","add_time","w_time","print_time",
+ "ocas_time","total_time"};
+ mwSize dims[2] = {1,1};
+
+ plhs[2] = mxCreateStructArray(2, dims, (sizeof(field_names)/sizeof(*field_names)),
+ field_names);
+
+ mxSetField(plhs[2],0,"nIter",mxCreateDoubleScalar((double)ocas.nIter));
+ mxSetField(plhs[2],0,"nCutPlanes",mxCreateDoubleScalar((double)ocas.nCutPlanes));
+ mxSetField(plhs[2],0,"nTrnErrors",mxCreateDoubleScalar(ocas.trn_err));
+ mxSetField(plhs[2],0,"Q_P",mxCreateDoubleScalar(ocas.Q_P));
+ mxSetField(plhs[2],0,"Q_D",mxCreateDoubleScalar(ocas.Q_D));
+ mxSetField(plhs[2],0,"init_time",mxCreateDoubleScalar(init_time));
+ mxSetField(plhs[2],0,"output_time",mxCreateDoubleScalar(ocas.output_time));
+ mxSetField(plhs[2],0,"sort_time",mxCreateDoubleScalar(ocas.sort_time));
+ mxSetField(plhs[2],0,"qp_solver_time",mxCreateDoubleScalar(ocas.qp_solver_time));
+ mxSetField(plhs[2],0,"add_time",mxCreateDoubleScalar(ocas.add_time));
+ mxSetField(plhs[2],0,"w_time",mxCreateDoubleScalar(ocas.w_time));
+ mxSetField(plhs[2],0,"print_time",mxCreateDoubleScalar(ocas.print_time));
+ mxSetField(plhs[2],0,"ocas_time",mxCreateDoubleScalar(ocas.ocas_time));
+ mxSetField(plhs[2],0,"total_time",mxCreateDoubleScalar(total_time));
+ mxSetField(plhs[2],0,"exitflag",mxCreateDoubleScalar((double)ocas.exitflag));
+
+ return;
+}
+
diff --git a/version.h b/version.h
index c79cd4a..5840f3c 100644
--- a/version.h
+++ b/version.h
@@ -1,3 +1,3 @@
/* this text will be printed in all programs */
-#define OCAS_VERSION "v0.95, (C) 2008-2012 Vojtech Franc, Soeren Sonnenburg"
+#define OCAS_VERSION "v0.97, 2013-02-26 (C) 2008-2012 Vojtech Franc, Soeren Sonnenburg"
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/libocas.git
More information about the debian-science-commits
mailing list