Permalink
Browse files

python wrapper around inner loop of partial_fit

  • Loading branch information...
1 parent a993319 commit 8e316c3916456062701294120c674729ab39d57a Claudio A Andreoni committed Dec 3, 2011
Showing with 352 additions and 253 deletions.
  1. +24 −23 asgd/Makefile
  2. +38 −69 asgd/asgd.c
  3. +3 −3 asgd/asgd.h
  4. +132 −0 asgd/asgd.py
  5. +11 −0 asgd/asgd_blas.h
  6. +103 −0 asgd/asgd_core.c
  7. +41 −0 asgd/asgd_core.h
  8. +0 −158 asgd/pasgd.pyx
View
@@ -1,39 +1,40 @@
CC = gcc
DEBUG = -g
-CFLAGS = -Wall -std=c99
+CFLAGS = -Wall -std=c99 -fPIC
-LIB = -lm -lcblas
INCDIR = -I/usr/include/
LIBDIR = -L/usr/lib64/atlas
-CYTHON = cython
-PFLAGS = $(CFLAGS) -pthread -fPIC -fwrapv -fno-strict-aliasing
-PINCDIR = $(INCDIR) -I/usr/include/python2.7 -I/usr/lib64/python2.7/site-packages/numpy/core/include
-PLIBDIR = $(LIBDIR)
-
-# decide which library to use for the linear algebra
-ifeq ($(TYPE),sse)
-LTYPE = -DASGD_SSE
-endif
ifeq ($(TYPE),blas)
-LTYPE = -DASGD_BLAS
-endif
-ifeq ($(TYPE),)
-LTYPE = simple_blas.c
+LIBS = -lm -lcblas
+DEFS = -DASGD_BLAS
+OBJS =
+else
+LIBS = -lm
+DEFS =
+OBJS = simple_blas.o
endif
-python_lib: pasgd.pyx asgd.c simple_blas.c
- $(CYTHON) pasgd.pyx
- $(CC) -shared $(PFLAGS) $(DEBUG) $(INCDIR) $(PINCDIR) $(LIBDIR) $(LIB) -o pasgd.so $(LTYPE) pasgd.c asgd.c
+asgd_unit: asgd.o asgd_core.o simple_blas.o tests/asgd_unit.c
+ $(CC) $(CFLAGS) $(DEBUG) $(INCDIR) $(LIBDIR) $(LIBS) $(DEFS) -o bin/asgd_unit $(OBJS) asgd.o asgd_core.o tests/asgd_unit.c
+
+asgd.o: simple_blas.o asgd_core.o asgd.c asgd.h
+ $(CC) $(CFLAGS) $(DEBUG) $(INCDIR) $(LIBDIR) $(LIBS) $(DEFS) -c -o asgd.o asgd.c
+
+asgd_core.o: simple_blas.o asgd_core.c asgd_core.h
+ $(CC) $(CFLAGS) $(DEBUG) $(INCDIR) $(LIBDIR) $(LIBS) $(DEFS) -c -o asgd_core.o asgd_core.c
+ $(CC) -shared $(CFLAGS) $(DEBUG) $(INCDIR) $(LIBDIR) $(LIBS) $(DEFS) -o asgd_core.so $(OBJS) asgd_core.c
-asgd_unit: tests/asgd_unit.c asgd.c simple_blas.c
- $(CC) $(CFLAGS) $(DEBUG) $(INCDIR) $(LIBDIR) $(LIB) -o bin/asgd_unit $(LTYPE) tests/asgd_unit.c asgd.c
+simple_blas_unit: simple_blas.c simple_blas.h tests/simple_blas_unit.c
+ $(CC) $(CFLAGS) $(DEBUG) $(INCDIR) $(LIBDIR) $(LIBS) $(DEFS) -o bin/simple_blas_unit tests/simple_blas_unit.c
-simple_blas_unit: tests/simple_blas_unit.c simple_blas.c
- $(CC) $(CFLAGS) $(DEBUG) $(INCDIR) $(LIBDIR) $(LIB) -o bin/simple_blas_unit tests/simple_blas_unit.c simple_blas.c
+simple_blas.o: simple_blas.c simple_blas.h
+ $(CC) $(CFLAGS) $(DEBUG) $(INCDIR) $(LIBDIR) $(LIBS) $(DEFS) -c -o simple_blas.o simple_blas.c
-# clean up the compilation byproducts
+# clean up compilation byproducts
.PHONY: clean
clean:
rm -fR *.o
+ rm -fR *.so
rm -f bin/asgd_unit
+ rm -f bin/simple_blas_unit
View
@@ -1,12 +1,7 @@
#include "asgd.h"
-#if defined ASGD_BLAS
-#include <cblas.h>
-#elif defined ASGD_SSE
-#include "sse_blas.h"
-#else
-#include "simple_blas.h"
-#endif
+#include "asgd_blas.h"
+#include "asgd_core.h"
matrix_t *matrix_init(
size_t rows,
@@ -108,10 +103,10 @@ void matrix_row_shuffle(matrix_t *m, int *r)
* Constructor for the Binary ASGD structure
*/
nb_asgd_t *nb_asgd_init(
- uint64_t n_feats,
+ long n_feats,
float sgd_step_size0,
float l2_reg,
- uint64_t n_iters,
+ long n_iters,
bool feedback)
{
nb_asgd_t *data = malloc(sizeof(*data));
@@ -158,65 +153,39 @@ void partial_fit(
matrix_t *X,
matrix_t *y)
{
-
- for (size_t i = 0; i < X->rows; ++i) {
-
- // compute margin //
- // TODO sgd_weights will become a matrix
- // notice that each row in X is also a column because of the stride
- float margin = matrix_get(y, i, 0) *
- cblas_sdsdot(
- X->cols,
- matrix_get(data->sgd_bias, 0, 0),
- matrix_row(X, i), 1,
- data->sgd_weights->data, 1);
-
- // update sgd //
- if (data->l2_reg != 0)
- {
- // TODO sgd_weights will become a matrix
- cblas_sscal(data->sgd_weights->rows,
- 1 - data->l2_reg * data->sgd_step_size,
- data->sgd_weights->data, 1);
- }
-
- if (margin < 1)
- {
- // TODO sgd_weights will become a matrix
- // TODO may be faster to leave sgd_weights on the stack
- cblas_saxpy(
- data->sgd_weights->rows,
- data->sgd_step_size * matrix_get(y, i, 0),
- matrix_row(X, i), 1,
- data->sgd_weights->data, 1);
-
- // TODO sgd_bias will become a vector
- matrix_set(data->sgd_bias, 0, 0,
- data->sgd_step_size * matrix_get(y, i, 0));
- }
-
- // update asgd //
- //matrix_t *asgd_weights = matrix_clone(data->asgd_weights);
- cblas_sscal(data->asgd_weights->rows,
- 1 - data->asgd_step_size,
- data->asgd_weights->data, 1);
- cblas_saxpy(data->asgd_weights->rows,
- data->asgd_step_size,
- data->sgd_weights->data, 1,
- data->asgd_weights->data, 1);
-
- matrix_set(data->asgd_bias, 0, 0,
- (1 - data->asgd_step_size) * matrix_get(data->asgd_bias, 0, 0) +
- data->asgd_step_size * matrix_get(data->sgd_bias, 0, 0));
-
- // update step_sizes //
- data->n_observs += 1;
- float sgd_step_size_scheduling = 1 + data->sgd_step_size0 * data->n_observs
- * data->sgd_step_size_scheduling_mul;
- data->sgd_step_size = data->sgd_step_size0 /
- pow(sgd_step_size_scheduling, data->sgd_step_size_scheduling_exp);
- data->asgd_step_size = 1.0f / data->n_observs;
- }
+ core_partial_fit(
+ &data->n_observs,
+ &data->sgd_step_size,
+ &data->asgd_step_size,
+
+ data->l2_reg,
+ data->sgd_step_size0,
+ data->sgd_step_size_scheduling_exp,
+ data->sgd_step_size_scheduling_mul,
+
+ data->sgd_weights->data,
+ data->sgd_weights->rows,
+ data->sgd_weights->cols,
+
+ data->sgd_bias->data,
+ data->sgd_bias->rows,
+ data->sgd_bias->cols,
+
+ data->asgd_weights->data,
+ data->asgd_weights->rows,
+ data->asgd_weights->cols,
+
+ data->asgd_bias->data,
+ data->asgd_bias->rows,
+ data->asgd_bias->cols,
+
+ X->data,
+ X->rows,
+ X->cols,
+
+ y->data,
+ y->rows,
+ y->cols);
}
void fit(
@@ -228,7 +197,7 @@ void fit(
mex_assert(X->rows > 1, "fit: X should be a matrix");
mex_assert(y->cols == 1, "fit: y should be a column vector");
- for (uint64_t i = 0; i < data->n_iters; ++i)
+ for (long i = 0; i < data->n_iters; ++i)
{
matrix_t *Xb = matrix_clone(X);
matrix_row_shuffle(Xb, r+i*Xb->rows);
View
@@ -57,14 +57,14 @@ struct nb_asgd
float asgd_step_size;
float asgd_step_size0;
- uint64_t n_observs;
+ long n_observs;
};
nb_asgd_t *nb_asgd_init(
- uint64_t n_feats,
+ long n_feats,
float sgd_step_size0,
float l2_reg,
- uint64_t n_iters,
+ long n_iters,
bool feedback);
void nb_asgd_destr(
View
@@ -0,0 +1,132 @@
+"""Averaging Stochastic Gradient Descent Classifier
+
+naive, non-optimized implementation
+"""
+import ctypes as ct
+import numpy as np
+from numpy import dot
+from itertools import izip
+
+
+class ASGD(object):
+
+ def __init__(self, n_features, sgd_step_size0=1e-2, l2_regularization=1e-3,
+ n_iterations=10, feedback=False, dtype=np.float32):
+
+ self.n_features = n_features
+ self.n_iterations = n_iterations
+ self.feedback = feedback
+
+ assert l2_regularization > 0
+ self.l2_regularization = l2_regularization
+ self.dtype = dtype
+
+ self.sgd_weights = np.zeros((n_features), dtype=dtype)
+ self.sgd_bias = np.zeros((1), dtype=dtype)
+ self.sgd_step_size0 = sgd_step_size0
+ self.sgd_step_size = sgd_step_size0
+ self.sgd_step_size_scheduling_exponent = 2. / 3
+ self.sgd_step_size_scheduling_multiplier = l2_regularization
+
+ self.asgd_weights = np.zeros((n_features), dtype=dtype)
+ self.asgd_bias = np.zeros((1), dtype=dtype)
+ self.asgd_step_size0 = 1
+ self.asgd_step_size = self.asgd_step_size0
+
+ self.n_observations = 0
+ self.core_lib = ct.CDLL("./asgd_core.so")
+
+
+ def partial_fit(self, X, y):
+
+ input_req = ['A', 'O', 'W', 'C']
+ output_req = ['A', 'O', 'W', 'C']
+ np.require(self.sgd_weights, dtype=np.float32, requirements=output_req)
+ np.require(self.sgd_bias, dtype=np.float32, requirements=output_req)
+ np.require(self.asgd_weights, dtype=np.float32, requirements=output_req)
+ np.require(self.asgd_bias, dtype=np.float32, requirements=output_req)
+ np.require(X, dtype=np.float32, requirements=input_req)
+ np.require(y, dtype=np.float32, requirements=input_req)
+
+ sgd_step_size0 = ct.c_float(self.sgd_step_size0)
+ sgd_step_size = ct.c_float(self.sgd_step_size)
+ sgd_step_size_scheduling_exponent = \
+ ct.c_float(self.sgd_step_size_scheduling_exponent)
+ sgd_step_size_scheduling_multiplier = \
+ ct.c_float(self.sgd_step_size_scheduling_multiplier)
+ sgd_weights = self.sgd_weights
+ sgd_bias = self.sgd_bias
+
+ asgd_weights = self.asgd_weights
+ asgd_bias = self.asgd_bias
+ asgd_step_size = ct.c_float(self.asgd_step_size)
+
+ l2_regularization = ct.c_float(self.l2_regularization)
+
+ n_observations = ct.c_long(self.n_observations)
+
+ self.core_lib.core_partial_fit(
+ ct.byref(n_observations),\
+ ct.byref(sgd_step_size),\
+ ct.byref(asgd_step_size),\
+ l2_regularization,\
+ sgd_step_size0,\
+ sgd_step_size_scheduling_exponent,\
+ sgd_step_size_scheduling_multiplier,\
+ sgd_weights.ctypes.data_as(ct.POINTER(ct.c_float)),\
+ ct.c_size_t(sgd_weights.shape[0]),\
+ ct.c_size_t(1),\
+ sgd_bias.ctypes.data_as(ct.POINTER(ct.c_float)),\
+ ct.c_size_t(1),\
+ ct.c_size_t(1),\
+ asgd_weights.ctypes.data_as(ct.POINTER(ct.c_float)),\
+ ct.c_size_t(asgd_weights.shape[0]),\
+ ct.c_size_t(1),\
+ asgd_bias.ctypes.data_as(ct.POINTER(ct.c_float)),\
+ ct.c_size_t(1),\
+ ct.c_size_t(1),\
+ X.ctypes.data_as(ct.POINTER(ct.c_float)),\
+ ct.c_size_t(X.shape[0]),\
+ ct.c_size_t(X.shape[1]),\
+ y.ctypes.data_as(ct.POINTER(ct.c_float)),\
+ ct.c_size_t(y.shape[0]),\
+ ct.c_size_t(y.shape[1]))
+
+ # --
+ self.sgd_weights = sgd_weights
+ self.sgd_bias = sgd_bias
+ self.sgd_step_size = sgd_step_size.value
+
+ self.asgd_weights = asgd_weights
+ self.asgd_bias = asgd_bias
+ self.asgd_step_size = asgd_step_size.value
+
+ self.n_observations = n_observations.value
+
+ def fit(self, X, y):
+
+ assert X.ndim == 2
+ assert y.ndim == 1
+
+ n_points, n_features = X.shape
+ assert n_features == self.n_features
+ assert n_points == y.size
+
+ n_iterations = self.n_iterations
+
+ for i in xrange(n_iterations):
+
+ idx = np.random.permutation(n_points)
+ Xb = X[idx]
+ yb = y[idx]
+ self.partial_fit(Xb, yb)
+
+ if self.feedback:
+ self.sgd_weights = self.asgd_weights
+ self.sgd_bias = self.asgd_bias
+
+ def decision_function(self, X):
+ return dot(self.asgd_weights, X.T) + self.asgd_bias
+
+ def predict(self, X):
+ return np.sign(self.decision_function(X))
View
@@ -0,0 +1,11 @@
+
+// determine which BLAS implementation to use
+
+#if defined ASGD_BLAS
+#include <cblas.h>
+#elif defined ASGD_SSE
+#include "sse_blas.h"
+#else
+#include "simple_blas.h"
+#endif
+
Oops, something went wrong.

0 comments on commit 8e316c3

Please sign in to comment.