|
Skylark (Sketching Library)
0.1
|
00001 /* 00002 * prox_cross_entropy.c 00003 * 00004 * Created on: Dec 10, 2013 00005 * Author: vikas 00006 */ 00007 00008 #include "Python.h" 00009 #include "arrayobject.h" 00010 #include "libproxoperators.h" 00011 #include <math.h> 00012 00013 static PyMethodDef _libproxoperatorsMethods[] = { 00014 {"crossentropy_prox", crossentropy_prox, METH_VARARGS}, 00015 {"crossentropy_obj", crossentropy_obj, METH_VARARGS}, 00016 {"hinge_prox", hinge_prox, METH_VARARGS}, 00017 {"hinge_obj", hinge_obj, METH_VARARGS}, 00018 {NULL, NULL} /* Sentinel - marks the end of this structure */ 00019 }; 00020 00021 void init_libproxoperators() { 00022 (void) Py_InitModule("_libproxoperators", _libproxoperatorsMethods); 00023 import_array(); // Must be present for NumPy. Called first after above line. 00024 } 00025 00026 static PyObject *crossentropy_prox(PyObject *self, PyObject *args) 00027 { 00028 PyArrayObject *x, *v, *y; 00029 double lambda, epsilon; 00030 int MAXITER, m, n, flag, i, DISPLAY; // ncomps=n*m=total number of matrix components in mat1 00031 00032 /* Parse tuples separately since args will differ between C fcns */ 00033 if (!PyArg_ParseTuple(args, "OOdOidi", 00034 &y, &v, &lambda, &x, &MAXITER, &epsilon,&DISPLAY)) return NULL; 00035 if (NULL == v) return NULL; 00036 00037 /* Check that object input is 'double' type and a matrix 00038 Not needed if python wrapper function checks before call to this routine */ 00039 if (v->descr->type_num != NPY_DOUBLE || v->nd != 2) { 00040 PyErr_SetString(PyExc_ValueError, 00041 "In not_doublematrix: array must be of type Float and 2 dimensional (m x n)."); 00042 return NULL; } 00043 00044 // printf("y: %d, %d", y->descr->type_num, y->nd); 00045 // if (y->descr->type_num != NPY_DOUBLE || y->nd != 2) { 00046 // PyErr_SetString(PyExc_ValueError, 00047 // "In not_doublematrix: array must be of type Float and 2 dimensional (m x 1)."); 00048 // return NULL; } 00049 00050 /* Get the dimensions of the input */ 00051 00052 PyObject *v_array = PyArray_FROM_OTF((PyObject *) v, NPY_DOUBLE, NPY_IN_ARRAY); 00053 PyObject *y_array = PyArray_FROM_OTF((PyObject *) y, NPY_DOUBLE, NPY_IN_ARRAY); 00054 PyObject *x_array = PyArray_FROM_OTF((PyObject *) x, NPY_DOUBLE, NPY_INOUT_ARRAY); 00055 00056 if (x_array == NULL || v_array == NULL) { 00057 Py_XDECREF(x_array); 00058 Py_XDECREF(y_array); 00059 Py_XDECREF(v_array); 00060 return NULL; 00061 } 00062 00063 // printf("Number of dimensions=%d\n", PyArray_NDIM(v_array)); 00064 00065 m = (int) PyArray_DIM(v_array,0); 00066 n = (int) PyArray_DIM(v_array,1); 00067 00068 double *cx = (double*)PyArray_DATA(x_array); 00069 double *cy = (double*) PyArray_DATA(y_array); 00070 double *cv = (double*)PyArray_DATA(v_array); 00071 00072 // printf("m=%d, n=%d,lambda=%f,MAXITER=%d,epsilon=%f\n", m, n, lambda, MAXITER, epsilon); 00073 // for(i=0;i<m;i++) 00074 // printf("%f ", cy[i]); 00075 // printf("\n"); 00076 00077 for(i=0;i<m;i++) 00078 flag = logexp((int) cy[i], cv + i*n, n, lambda, cx + i*n, MAXITER, epsilon, DISPLAY); 00079 00080 Py_DECREF(x_array); 00081 Py_DECREF(v_array); 00082 Py_DECREF(y_array); 00083 00084 // printf("returned %d", flag); 00085 00086 PyObject *ret = Py_BuildValue("i", flag); 00087 return ret; 00088 } 00089 00090 00091 static PyObject *hinge_prox(PyObject *self, PyObject *args) 00092 { 00093 PyArrayObject *x, *v, *y; 00094 double lambda, yv, yy, *cvv, *cxx; 00095 int m, n, i, j,label; // ncomps=n*m=total number of matrix components in mat1 00096 00097 /* Parse tuples separately since args will differ between C fcns */ 00098 if (!PyArg_ParseTuple(args, "OOdO", 00099 &y, &v, &lambda, &x)) return NULL; 00100 if (NULL == v) return NULL; 00101 00102 /* Check that object input is 'double' type and a matrix 00103 Not needed if python wrapper function checks before call to this routine */ 00104 if (v->descr->type_num != NPY_DOUBLE || v->nd != 2) { 00105 PyErr_SetString(PyExc_ValueError, 00106 "In not_doublematrix: array must be of type Float and 2 dimensional (m x n)."); 00107 return NULL; } 00108 00109 // printf("y: %d, %d", y->descr->type_num, y->nd); 00110 // if (y->descr->type_num != NPY_DOUBLE || y->nd != 2) { 00111 // PyErr_SetString(PyExc_ValueError, 00112 // "In not_doublematrix: array must be of type Float and 2 dimensional (m x 1)."); 00113 // return NULL; } 00114 00115 /* Get the dimensions of the input */ 00116 00117 PyObject *v_array = PyArray_FROM_OTF((PyObject *) v, NPY_DOUBLE, NPY_IN_ARRAY); 00118 PyObject *y_array = PyArray_FROM_OTF((PyObject *) y, NPY_DOUBLE, NPY_IN_ARRAY); 00119 PyObject *x_array = PyArray_FROM_OTF((PyObject *) x, NPY_DOUBLE, NPY_INOUT_ARRAY); 00120 00121 if (x_array == NULL || v_array == NULL) { 00122 Py_XDECREF(x_array); 00123 Py_XDECREF(y_array); 00124 Py_XDECREF(v_array); 00125 return NULL; 00126 } 00127 00128 // printf("Number of dimensions=%d\n", PyArray_NDIM(v_array)); 00129 00130 m = (int) PyArray_DIM(v_array,0); 00131 n = (int) PyArray_DIM(v_array,1); 00132 00133 double *cx = (double*)PyArray_DATA(x_array); 00134 double *cy = (double*) PyArray_DATA(y_array); 00135 double *cv = (double*)PyArray_DATA(v_array); 00136 00137 // printf("m=%d, n=%d,lambda=%f,MAXITER=%d,epsilon=%f\n", m, n, lambda, MAXITER, epsilon); 00138 // for(i=0;i<m;i++) 00139 // printf("%f ", cy[i]); 00140 // printf("\n"); 00141 //flag = logexp((int) cy[i], cv + i*n, n, lambda, cx + i*n); 00142 if(n==1) { // We assume cy has +1 or -1 entries for n=1 outputs 00143 for(i=0;i<m;i++) { 00144 yv = cy[i]*cv[i]; 00145 if (yv>1.0) 00146 cx[i] = cv[i]; 00147 else { 00148 if(yv<1.0-lambda) 00149 cx[i] = cv[i] + lambda*cy[i]; 00150 else 00151 cx[i] = cy[i]; 00152 } 00153 } 00154 } 00155 00156 if (n>1) { 00157 for(i=0;i<m;i++) { 00158 label = (int) cy[i]; 00159 cvv = cv + i*n; 00160 cxx = cx + i*n; 00161 for(j=0;j<n;j++) { 00162 yv = cvv[j]; 00163 yy = +1.0; 00164 if(!(j==label)) { 00165 yv = -yv; 00166 yy = -1.0; 00167 } 00168 if (yv>1.0) 00169 cxx[j] = cvv[j]; 00170 else { 00171 if(yv<1.0-lambda) 00172 cxx[j] = cvv[j] + lambda*yy; 00173 else 00174 cxx[j] = yy; 00175 } 00176 } 00177 } 00178 } 00179 00180 Py_DECREF(x_array); 00181 Py_DECREF(v_array); 00182 Py_DECREF(y_array); 00183 00184 // printf("returned %d", flag); 00185 00186 PyObject *ret = Py_BuildValue("i", 1.0); 00187 return ret; 00188 } 00189 00190 00191 int logexp(int index, double* v, int n, double lambda, double* x, int MAXITER, double epsilon, int DISPLAY) { 00192 /* solution to - log exp(x(i))/sum(exp(x(j))) + lambda/2 ||x - v||_2^2 */ 00193 /* n is length of v and x */ 00194 /* writes over x */ 00195 double alpha = 0.1; 00196 double beta = 0.5; 00197 int iter, i; 00198 double t, logsum, p, pu, pptil, decrement; 00199 double *u = (double *) malloc(n*sizeof(double)); 00200 double *z = (double *) malloc(n*sizeof(double)); 00201 double *grad = (double *) malloc(n*sizeof(double)); 00202 double newobj=0.0, obj=0.0; 00203 obj = objective(index, x, v, n, lambda); 00204 00205 for(iter=0;iter<MAXITER;iter++) { 00206 logsum = logsumexp(x,n); 00207 if(DISPLAY) 00208 printf("iter=%d, obj=%f\n", iter, obj); 00209 pu = 0.0; 00210 pptil = 0.0; 00211 for(i=0;i<n;i++) { 00212 p = exp(x[i] - logsum); 00213 grad[i] = p + lambda*(x[i] - v[i]); 00214 if(i==index) 00215 grad[i] += -1.0; 00216 u[i] = grad[i]/(p+lambda); 00217 pu += p*u[i]; 00218 z[i] = p/(p+lambda); 00219 pptil += z[i]*p; 00220 } 00221 decrement = 0.0; 00222 for(i=0;i<n;i++) { 00223 u[i] -= (pu/pptil)*z[i]; 00224 decrement += grad[i]*u[i]; 00225 } 00226 if (decrement < 2*epsilon) { 00227 free(u); 00228 free(z); 00229 free(grad); 00230 return 0; 00231 } 00232 t = 1.0; 00233 while(1) { 00234 for(i=0;i<n;i++) 00235 z[i] = x[i] - t*u[i]; 00236 newobj = objective(index, z, v, n, lambda); 00237 if (newobj <= obj + alpha*t*decrement) 00238 break; 00239 t = beta*t; 00240 } 00241 for(i=0;i<n;i++) 00242 x[i] = z[i]; 00243 obj = newobj; 00244 } 00245 free(u); 00246 free(z); 00247 free(grad); 00248 return 1; 00249 } 00250 00251 00252 double objective(int index, double* x, double* v, int n, double lambda) { 00253 double nrmsqr = normsquare(x,v,n); 00254 double obj = -x[index] + logsumexp(x, n) + 0.5*lambda*nrmsqr; 00255 return obj; 00256 } 00257 00258 double normsquare(double* x, double* y, int n){ 00259 double nrm = 0.0; 00260 int i; 00261 for(i=0;i<n;i++) 00262 nrm+= pow(x[i] - y[i], 2); 00263 return nrm; 00264 } 00265 00266 double logsumexp(double* x, int n) { 00267 int i; 00268 double max=x[0]; 00269 double f = 0.0; 00270 for(i=0;i<n;i++) { 00271 if (x[i]>max) { 00272 max = x[i]; 00273 } 00274 } 00275 for(i=0;i<n;i++) 00276 f += exp(x[i] - max); 00277 f = max + log(f); 00278 00279 return f; 00280 } 00281 00282 00283 static PyObject *crossentropy_obj(PyObject *self, PyObject *args) 00284 { 00285 PyArrayObject *x, *y; 00286 int m,n, i; 00287 double obj; 00288 /* Parse tuples separately since args will differ between C fcns */ 00289 if (!PyArg_ParseTuple(args, "OO", 00290 &y, &x)) return NULL; 00291 if (NULL == x) return NULL; 00292 00293 /* Check that object input is 'double' type and a matrix 00294 Not needed if python wrapper function checks before call to this routine */ 00295 if (x->descr->type_num != NPY_DOUBLE || x->nd != 2) { 00296 PyErr_SetString(PyExc_ValueError, 00297 "In not_doublematrix: array must be of type Float and 2 dimensional (m x n)."); 00298 return NULL; } 00299 00300 if (y->descr->type_num != NPY_DOUBLE || y->nd != 2) { 00301 PyErr_SetString(PyExc_ValueError, 00302 "In not_doublematrix: array must be of type Float and 2 dimensional (m x n)."); 00303 return NULL; } 00304 00305 // printf("y: %d, %d", y->descr->type_num, y->nd); 00306 // if (y->descr->type_num != NPY_DOUBLE || y->nd != 2) { 00307 // PyErr_SetString(PyExc_ValueError, 00308 // "In not_doublematrix: array must be of type Float and 2 dimensional (m x 1)."); 00309 // return NULL; } 00310 00311 /* Get the dimensions of the input */ 00312 00313 PyObject *y_array = PyArray_FROM_OTF((PyObject *) y, NPY_DOUBLE, NPY_IN_ARRAY); 00314 PyObject *x_array = PyArray_FROM_OTF((PyObject *) x, NPY_DOUBLE, NPY_IN_ARRAY); 00315 00316 if (x_array == NULL || y_array == NULL) { 00317 Py_XDECREF(x_array); 00318 Py_XDECREF(y_array); 00319 return NULL; 00320 } 00321 00322 // printf("Number of dimensions=%d\n", PyArray_NDIM(v_array)); 00323 00324 m = (int) PyArray_DIM(x_array,0); 00325 n = (int) PyArray_DIM(x_array,1); 00326 00327 double *cx = (double*)PyArray_DATA(x_array); 00328 double *cy = (double*) PyArray_DATA(y_array); 00329 00330 obj = 0.0; 00331 for(i=0;i<m;i++) 00332 obj += -*(cx + i*n + (int) cy[i]) + logsumexp(cx + i*n, n); 00333 00334 00335 Py_DECREF(x_array); 00336 Py_DECREF(y_array); 00337 00338 // printf("returned %d", flag); 00339 00340 PyObject *ret = Py_BuildValue("d", obj); 00341 return ret; 00342 } 00343 00344 static PyObject *hinge_obj(PyObject *self, PyObject *args) 00345 { 00346 PyArrayObject *x, *y; 00347 int m,n, i,j, label; 00348 double obj, yx, yy, *cxx; 00349 /* Parse tuples separately since args will differ between C fcns */ 00350 if (!PyArg_ParseTuple(args, "OO", 00351 &y, &x)) return NULL; 00352 if (NULL == x) return NULL; 00353 00354 /* Check that object input is 'double' type and a matrix 00355 Not needed if python wrapper function checks before call to this routine */ 00356 if (x->descr->type_num != NPY_DOUBLE || x->nd != 2) { 00357 PyErr_SetString(PyExc_ValueError, 00358 "In not_doublematrix: array must be of type Float and 2 dimensional (m x n)."); 00359 return NULL; } 00360 00361 if (y->descr->type_num != NPY_DOUBLE || y->nd != 2) { 00362 PyErr_SetString(PyExc_ValueError, 00363 "In not_doublematrix: array must be of type Float and 2 dimensional (m x n)."); 00364 return NULL; } 00365 00366 // printf("y: %d, %d", y->descr->type_num, y->nd); 00367 // if (y->descr->type_num != NPY_DOUBLE || y->nd != 2) { 00368 // PyErr_SetString(PyExc_ValueError, 00369 // "In not_doublematrix: array must be of type Float and 2 dimensional (m x 1)."); 00370 // return NULL; } 00371 00372 /* Get the dimensions of the input */ 00373 00374 PyObject *y_array = PyArray_FROM_OTF((PyObject *) y, NPY_DOUBLE, NPY_IN_ARRAY); 00375 PyObject *x_array = PyArray_FROM_OTF((PyObject *) x, NPY_DOUBLE, NPY_IN_ARRAY); 00376 00377 if (x_array == NULL || y_array == NULL) { 00378 Py_XDECREF(x_array); 00379 Py_XDECREF(y_array); 00380 return NULL; 00381 } 00382 00383 // printf("Number of dimensions=%d\n", PyArray_NDIM(v_array)); 00384 00385 m = (int) PyArray_DIM(x_array,0); 00386 n = (int) PyArray_DIM(x_array,1); 00387 00388 double *cx = (double*)PyArray_DATA(x_array); 00389 double *cy = (double*) PyArray_DATA(y_array); 00390 00391 obj = 0.0; 00392 00393 if(n==1) { 00394 for(i=0;i<m;i++) { 00395 yx = cx[i]*cy[i]; 00396 if(yx<1.0) 00397 obj += (1.0 - yx); 00398 } 00399 } 00400 00401 00402 if(n>1) { 00403 for(i=0;i<m;i++) { 00404 label = (int) cy[i]; 00405 cxx = cx + i*n; 00406 for(j=0;j<n;j++) { 00407 yy = -1.0; 00408 if (j==label) 00409 yy = +1.0; 00410 yx = cxx[j]*yy; 00411 if(yx<1.0) 00412 obj += (1.0 - yx); 00413 } 00414 } 00415 } 00416 00417 00418 Py_DECREF(x_array); 00419 Py_DECREF(y_array); 00420 00421 // printf("returned %d", flag); 00422 00423 PyObject *ret = Py_BuildValue("d", obj); 00424 return ret; 00425 }
1.7.6.1