pr_loqo.cpp

Go to the documentation of this file.
00001 /*
00002  * This program is free software; you can redistribute it and/or modify
00003  * it under the terms of the GNU General Public License as published by
00004  * the Free Software Foundation; either version 3 of the License, or
00005  * (at your option) any later version.
00006  *
00007  * Purpose:     solves quadratic programming problem for pattern recognition
00008  *              for support vectors
00009  *
00010  * Written (W) 1997-1998 Alex J. Smola
00011  * Written (W) 1999-2008 Soeren Sonnenburg
00012  * Written (W) 1999-2008 Gunnar Raetsch
00013  * Copyright (C) 1997-2008 Fraunhofer Institute FIRST and Max-Planck-Society
00014  */
00015 
00016 #include "lib/common.h"
00017 #include "lib/io.h"
00018 #include "lib/lapack.h"
00019 #include "lib/Mathematics.h"
00020 #include "classifier/svm/pr_loqo.h"
00021 
00022 #include <math.h>
00023 #include <time.h>
00024 #include <stdlib.h>
00025 #include <stdio.h>
00026 #ifdef HAVE_LAPACK
00027 extern "C" {
00028 #ifdef HAVE_ATLAS
00029 #include <clapack.h>
00030 #endif
00031 }
00032 #endif
00033 
00034 #define PREDICTOR 1
00035 #define CORRECTOR 2
00036 
00037 /*****************************************************************
00038   replace this by any other function that will exit gracefully
00039   in a larger system
00040   ***************************************************************/
00041 
00042 void nrerror(CHAR error_text[])
00043 {
00044     SG_SDEBUG("terminating optimizer - %s\n", error_text);
00045  // exit(1); 
00046 }
00047 
00048 /*****************************************************************
00049    taken from numerical recipes and modified to accept pointers
00050    moreover numerical recipes code seems to be buggy (at least the
00051    ones on the web)
00052 
00053    cholesky solver and backsubstitution
00054    leaves upper right triangle intact (rows first order)
00055    ***************************************************************/
00056 
00057 #ifdef HAVE_LAPACK
00058 bool choldc(double* a, int n, double* p)
00059 {
00060     if (n<=0)
00061         return false;
00062 
00063     double* a2=new double[n*n];
00064 
00065     for (int i=0; i<n; i++)
00066     {
00067         for (int j=0; j<n; j++)
00068             a2[n*i+j]=a[n*i+j];
00069     }
00070 
00071     int result=clapack_dpotrf(CblasRowMajor, CblasUpper, n, a2, n);
00072 
00073     for (int i=0; i<n; i++)
00074         p[i]=a2[(n+1)*i];
00075 
00076     for (int i=0; i<n; i++)
00077     {
00078         for (int j=i+1; j<n; j++)
00079         {
00080             a[n*j+i]=a2[n*i+j];
00081         }
00082     }
00083 
00084     if (result>0)
00085         SG_SDEBUG("Choldc failed, matrix not positive definite\n");
00086 
00087     delete[] a2;
00088     
00089     return result==0;
00090 }
00091 #else
00092 bool choldc(double a[], int n, double p[])
00093 {
00094     void nrerror(CHAR error_text[]);
00095     int i, j, k;
00096     double sum;
00097 
00098     for (i = 0; i < n; i++)
00099     {
00100         for (j = i; j < n; j++)
00101         {
00102             sum=a[n*i + j];
00103 
00104             for (k=i-1; k>=0; k--)
00105                 sum -= a[n*i + k]*a[n*j + k];
00106 
00107             if (i == j)
00108             {
00109                 if (sum <= 0.0)
00110                 {
00111                     SG_SERROR("Choldc failed, matrix not positive definite");
00112                     sum = 0.0;
00113                     return false;
00114                 }
00115 
00116                 p[i]=sqrt(sum);
00117 
00118             } 
00119             else
00120                 a[n*j + i] = sum/p[i];
00121         }
00122     }
00123 
00124     return true;
00125 }
00126 #endif
00127 
00128 void cholsb(double a[], int n, double p[], double b[], double x[])
00129 {
00130   int i, k;
00131   double sum;
00132 
00133   for (i=0; i<n; i++) {
00134     sum=b[i];
00135     for (k=i-1; k>=0; k--) sum -= a[n*i + k]*x[k];
00136     x[i]=sum/p[i];
00137   }
00138 
00139   for (i=n-1; i>=0; i--) {
00140     sum=x[i];
00141     for (k=i+1; k<n; k++) sum -= a[n*k + i]*x[k];
00142     x[i]=sum/p[i];
00143   }
00144 }
00145 
00146 /*****************************************************************
00147   sometimes we only need the forward or backward pass of the
00148   backsubstitution, hence we provide these two routines separately 
00149   ***************************************************************/
00150 
00151 void chol_forward(double a[], int n, double p[], double b[], double x[])
00152 {
00153   int i, k;
00154   double sum;
00155 
00156   for (i=0; i<n; i++) {
00157     sum=b[i];
00158     for (k=i-1; k>=0; k--) sum -= a[n*i + k]*x[k];
00159     x[i]=sum/p[i];
00160   }
00161 }
00162 
00163 void chol_backward(double a[], int n, double p[], double b[], double x[])
00164 {
00165   int i, k;
00166   double sum;
00167 
00168   for (i=n-1; i>=0; i--) {
00169     sum=b[i];
00170     for (k=i+1; k<n; k++) sum -= a[n*k + i]*x[k];
00171     x[i]=sum/p[i];
00172   }
00173 }
00174 
00175 /*****************************************************************
00176   solves the system | -H_x A' | |x_x| = |c_x|
00177                     |  A   H_y| |x_y|   |c_y|
00178 
00179   with H_x (and H_y) positive (semidefinite) matrices
00180   and n, m the respective sizes of H_x and H_y
00181 
00182   for variables see pg. 48 of notebook or do the calculations on a
00183   sheet of paper again
00184 
00185   predictor solves the whole thing, corrector assues that H_x didn't
00186   change and relies on the results of the predictor. therefore do
00187   _not_ modify workspace
00188 
00189   if you want to speed tune anything in the code here's the right
00190   place to do so: about 95% of the time is being spent in
00191   here. something like an iterative refinement would be nice,
00192   especially when switching from double to single precision. if you
00193   have a fast parallel cholesky use it instead of the numrec
00194   implementations.
00195 
00196   side effects: changes H_y (but this is just the unit matrix or zero anyway
00197   in our case)
00198   ***************************************************************/
00199 
00200 bool solve_reduced(int n, int m, double h_x[], double h_y[], 
00201            double a[], double x_x[], double x_y[],
00202            double c_x[], double c_y[],
00203            double workspace[], int step)
00204 {
00205   int i,j,k;
00206 
00207   double *p_x;
00208   double *p_y;
00209   double *t_a;
00210   double *t_c;
00211   double *t_y;
00212 
00213   p_x = workspace;      /* together n + m + n*m + n + m = n*(m+2)+2*m */
00214   p_y = p_x + n;
00215   t_a = p_y + m;
00216   t_c = t_a + n*m;
00217   t_y = t_c + n;
00218 
00219   if (step == PREDICTOR) {
00220     if (!choldc(h_x, n, p_x))   /* do cholesky decomposition */
00221         return false;
00222 
00223     for (i=0; i<m; i++)         /* forward pass for A' */
00224       chol_forward(h_x, n, p_x, a+i*n, t_a+i*n);
00225                 
00226     for (i=0; i<m; i++)         /* compute (h_y + a h_x^-1A') */
00227       for (j=i; j<m; j++)
00228     for (k=0; k<n; k++) 
00229       h_y[m*i + j] += t_a[n*j + k] * t_a[n*i + k];
00230                 
00231     choldc(h_y, m, p_y);    /* and cholesky decomposition */
00232   }
00233   
00234   chol_forward(h_x, n, p_x, c_x, t_c);
00235                 /* forward pass for c */
00236 
00237   for (i=0; i<m; i++) {     /* and solve for x_y */
00238     t_y[i] = c_y[i];
00239     for (j=0; j<n; j++)
00240       t_y[i] += t_a[i*n + j] * t_c[j];
00241   }
00242 
00243   cholsb(h_y, m, p_y, t_y, x_y);
00244 
00245   for (i=0; i<n; i++) {     /* finally solve for x_x */
00246     t_c[i] = -t_c[i];
00247     for (j=0; j<m; j++)
00248       t_c[i] += t_a[j*n + i] * x_y[j];
00249   }
00250 
00251   chol_backward(h_x, n, p_x, t_c, x_x);
00252   return true;
00253 }
00254 
00255 /*****************************************************************
00256   matrix vector multiplication (symmetric matrix but only one triangle
00257   given). computes m*x = y
00258   no need to tune it as it's only of O(n^2) but cholesky is of
00259   O(n^3). so don't waste your time _here_ although it isn't very
00260   elegant. 
00261   ***************************************************************/
00262 
00263 void matrix_vector(int n, double m[], double x[], double y[])
00264 {
00265   int i, j;
00266 
00267   for (i=0; i<n; i++) {
00268     y[i] = m[(n+1) * i] * x[i];
00269 
00270     for (j=0; j<i; j++)
00271       y[i] += m[i + n*j] * x[j];
00272 
00273     for (j=i+1; j<n; j++) 
00274       y[i] += m[n*i + j] * x[j]; 
00275   }
00276 }
00277 
00278 /*****************************************************************
00279   call only this routine; this is the only one you're interested in
00280   for doing quadratical optimization
00281 
00282   the restart feature exists but it may not be of much use due to the
00283   fact that an initial setting, although close but not very close the
00284   the actual solution will result in very good starting diagnostics
00285   (primal and dual feasibility and small infeasibility gap) but incur
00286   later stalling of the optimizer afterwards as we have to enforce
00287   positivity of the slacks.
00288   ***************************************************************/
00289 
00290 int pr_loqo(int n, int m, double c[], double h_x[], double a[], double b[],
00291         double l[], double u[], double primal[], double dual[], 
00292         int verb, double sigfig_max, int counter_max, 
00293         double margin, double bound, int restart) 
00294 {
00295   /* the knobs to be tuned ... */
00296   /* double margin = -0.95;    we will go up to 95% of the
00297                    distance between old variables and zero */
00298   /* double bound = 10;        preset value for the start. small
00299                    values give good initial
00300                    feasibility but may result in slow
00301                    convergence afterwards: we're too
00302                    close to zero */
00303   /* to be allocated */
00304   double *workspace;
00305   double *diag_h_x;
00306   double *h_y;
00307   double *c_x;
00308   double *c_y;
00309   double *h_dot_x;
00310   double *rho;
00311   double *nu;
00312   double *tau;
00313   double *sigma;
00314   double *gamma_z;
00315   double *gamma_s;  
00316 
00317   double *hat_nu;
00318   double *hat_tau;
00319 
00320   double *delta_x;
00321   double *delta_y;
00322   double *delta_s;
00323   double *delta_z;
00324   double *delta_g;
00325   double *delta_t;
00326 
00327   double *d;
00328 
00329   /* from the header - pointers into primal and dual */
00330   double *x;
00331   double *y;
00332   double *g;
00333   double *z;
00334   double *s;
00335   double *t;  
00336 
00337   /* auxiliary variables */
00338   double b_plus_1;
00339   double c_plus_1;
00340 
00341   double x_h_x;
00342   double primal_inf;
00343   double dual_inf;
00344 
00345   double sigfig;
00346   double primal_obj, dual_obj;
00347   double mu;
00348   double alfa=-1;
00349   int counter = 0;
00350 
00351   int status = STILL_RUNNING;
00352 
00353   int i,j;
00354 
00355   /* memory allocation */
00356   workspace = (double*) malloc((n*(m+2)+2*m)*sizeof(double));
00357   diag_h_x  = (double*) malloc(n*sizeof(double));
00358   h_y       = (double*) malloc(m*m*sizeof(double));
00359   c_x       = (double*) malloc(n*sizeof(double));
00360   c_y       = (double*) malloc(m*sizeof(double));
00361   h_dot_x   = (double*) malloc(n*sizeof(double));
00362 
00363   rho       = (double*) malloc(m*sizeof(double));
00364   nu        = (double*) malloc(n*sizeof(double));
00365   tau       = (double*) malloc(n*sizeof(double));
00366   sigma     = (double*) malloc(n*sizeof(double));
00367 
00368   gamma_z   = (double*) malloc(n*sizeof(double));
00369   gamma_s   = (double*) malloc(n*sizeof(double));
00370 
00371   hat_nu    = (double*) malloc(n*sizeof(double));
00372   hat_tau   = (double*) malloc(n*sizeof(double));
00373 
00374   delta_x   = (double*) malloc(n*sizeof(double));
00375   delta_y   = (double*) malloc(m*sizeof(double));
00376   delta_s   = (double*) malloc(n*sizeof(double));
00377   delta_z   = (double*) malloc(n*sizeof(double));
00378   delta_g   = (double*) malloc(n*sizeof(double));
00379   delta_t   = (double*) malloc(n*sizeof(double));
00380 
00381   d         = (double*) malloc(n*sizeof(double));
00382 
00383   /* pointers into the external variables */
00384   x = primal;           /* n */
00385   g = x + n;            /* n */
00386   t = g + n;            /* n */
00387 
00388   y = dual;         /* m */
00389   z = y + m;            /* n */
00390   s = z + n;            /* n */
00391 
00392   /* initial settings */
00393   b_plus_1 = 1;
00394   c_plus_1 = 0;
00395   for (i=0; i<n; i++) c_plus_1 += c[i];
00396 
00397   /* get diagonal terms */
00398   for (i=0; i<n; i++) diag_h_x[i] = h_x[(n+1)*i]; 
00399   
00400   /* starting point */
00401   if (restart == 1) {
00402                 /* x, y already preset */
00403     for (i=0; i<n; i++) {   /* compute g, t for primal feasibility */
00404       g[i] = CMath::max(CMath::abs(x[i] - l[i]), bound);
00405       t[i] = CMath::max(CMath::abs(u[i] - x[i]), bound); 
00406     }
00407 
00408     matrix_vector(n, h_x, x, h_dot_x); /* h_dot_x = h_x * x */
00409 
00410     for (i=0; i<n; i++) {   /* sigma is a dummy variable to calculate z, s */
00411       sigma[i] = c[i] + h_dot_x[i];
00412       for (j=0; j<m; j++)
00413     sigma[i] -= a[n*j + i] * y[j];
00414 
00415       if (sigma[i] > 0) {
00416     s[i] = bound;
00417     z[i] = sigma[i] + bound;
00418       }
00419       else {
00420     s[i] = bound - sigma[i];
00421     z[i] = bound;
00422       }
00423     }
00424   }
00425   else {            /* use default start settings */
00426     for (i=0; i<m; i++)
00427       for (j=i; j<m; j++)
00428     h_y[i*m + j] = (i==j) ? 1 : 0;
00429     
00430     for (i=0; i<n; i++) {
00431       c_x[i] = c[i];
00432       h_x[(n+1)*i] += 1;
00433     }
00434     
00435     for (i=0; i<m; i++)
00436       c_y[i] = b[i];
00437     
00438     /* and solve the system [-H_x A'; A H_y] [x, y] = [c_x; c_y] */
00439     solve_reduced(n, m, h_x, h_y, a, x, y, c_x, c_y, workspace,
00440           PREDICTOR);
00441     
00442     /* initialize the other variables */
00443     for (i=0; i<n; i++) {
00444       g[i] = CMath::max(CMath::abs(x[i] - l[i]), bound);
00445       z[i] = CMath::max(CMath::abs(x[i]), bound);
00446       t[i] = CMath::max(CMath::abs(u[i] - x[i]), bound); 
00447       s[i] = CMath::max(CMath::abs(x[i]), bound); 
00448     }
00449   }
00450 
00451   for (i=0, mu=0; i<n; i++)
00452     mu += z[i] * g[i] + s[i] * t[i];
00453   mu = mu / (2*n);
00454 
00455   /* the main loop */
00456   if (verb >= STATUS) {
00457       SG_SDEBUG("counter | pri_inf  | dual_inf  | pri_obj   | dual_obj  | ");
00458       SG_SDEBUG("sigfig | alpha  | nu \n");
00459       SG_SDEBUG("-------------------------------------------------------");
00460       SG_SDEBUG("---------------------------\n");
00461   }
00462   
00463   while (status == STILL_RUNNING) {
00464     /* predictor */
00465     
00466     /* put back original diagonal values */
00467     for (i=0; i<n; i++) 
00468       h_x[(n+1) * i] = diag_h_x[i];
00469 
00470     matrix_vector(n, h_x, x, h_dot_x); /* compute h_dot_x = h_x * x */
00471 
00472     for (i=0; i<m; i++) {
00473       rho[i] = b[i];
00474       for (j=0; j<n; j++)
00475     rho[i] -= a[n*i + j] * x[j];
00476     }
00477     
00478     for (i=0; i<n; i++) {
00479       nu[i] = l[i] - x[i] + g[i];
00480       tau[i] = u[i] - x[i] - t[i];
00481 
00482       sigma[i] = c[i] - z[i] + s[i] + h_dot_x[i];
00483       for (j=0; j<m; j++)
00484     sigma[i] -= a[n*j + i] * y[j];
00485 
00486       gamma_z[i] = - z[i];
00487       gamma_s[i] = - s[i];
00488     }
00489 
00490     /* instrumentation */
00491     x_h_x = 0;
00492     primal_inf = 0;
00493     dual_inf = 0;
00494     
00495     for (i=0; i<n; i++) {
00496       x_h_x += h_dot_x[i] * x[i];
00497       primal_inf += CMath::sq(tau[i]);
00498       primal_inf += CMath::sq(nu[i]);
00499       dual_inf += CMath::sq(sigma[i]);
00500     }
00501     for (i=0; i<m; i++) 
00502       primal_inf += CMath::sq(rho[i]);
00503     primal_inf = sqrt(primal_inf)/b_plus_1;
00504     dual_inf = sqrt(dual_inf)/c_plus_1;
00505 
00506     primal_obj = 0.5 * x_h_x;
00507     dual_obj = -0.5 * x_h_x;
00508     for (i=0; i<n; i++) {
00509       primal_obj += c[i] * x[i];
00510       dual_obj += l[i] * z[i] - u[i] * s[i];
00511     }
00512     for (i=0; i<m; i++)
00513       dual_obj += b[i] * y[i];
00514 
00515     sigfig = log10(CMath::abs(primal_obj) + 1) -
00516              log10(CMath::abs(primal_obj - dual_obj));
00517     sigfig = CMath::max(sigfig, 0.0);
00518            
00519     /* the diagnostics - after we computed our results we will
00520        analyze them */
00521 
00522     if (counter > counter_max) status = ITERATION_LIMIT;
00523     if (sigfig  > sigfig_max)  status = OPTIMAL_SOLUTION;
00524     if (primal_inf > 10e100)   status = PRIMAL_INFEASIBLE;
00525     if (dual_inf > 10e100)     status = DUAL_INFEASIBLE;
00526     if ((primal_inf > 10e100) & (dual_inf > 10e100)) status = PRIMAL_AND_DUAL_INFEASIBLE;
00527     if (CMath::abs(primal_obj) > 10e100) status = PRIMAL_UNBOUNDED;
00528     if (CMath::abs(dual_obj) > 10e100) status = DUAL_UNBOUNDED;
00529 
00530     /* write some nice routine to enforce the time limit if you
00531        _really_ want, however it's quite useless as you can compute
00532        the time from the maximum number of iterations as every
00533        iteration costs one cholesky decomposition plus a couple of 
00534        backsubstitutions */
00535 
00536     /* generate report */
00537     if ((verb >= FLOOD) | ((verb == STATUS) & (status != 0)))
00538      SG_SDEBUG("%7i | %.2e | %.2e | % .2e | % .2e | %6.3f | %.4f | %.2e\n",
00539          counter, primal_inf, dual_inf, primal_obj, dual_obj,
00540          sigfig, alfa, mu);
00541 
00542     counter++;
00543 
00544     if (status == 0) {      /* we may keep on going, otherwise
00545                    it'll cost one loop extra plus a
00546                    messed up main diagonal of h_x */
00547       /* intermediate variables (the ones with hat) */
00548       for (i=0; i<n; i++) {
00549     hat_nu[i] = nu[i] + g[i] * gamma_z[i] / z[i];
00550     hat_tau[i] = tau[i] - t[i] * gamma_s[i] / s[i];
00551     /* diagonal terms */
00552     d[i] = z[i] / g[i] + s[i] / t[i];
00553       }
00554       
00555       /* initialization before the cholesky solver */
00556       for (i=0; i<n; i++) {
00557     h_x[(n+1)*i] = diag_h_x[i] + d[i];
00558     c_x[i] = sigma[i] - z[i] * hat_nu[i] / g[i] - 
00559       s[i] * hat_tau[i] / t[i];
00560       }
00561       for (i=0; i<m; i++) {
00562     c_y[i] = rho[i];
00563     for (j=i; j<m; j++) 
00564       h_y[m*i + j] = 0;
00565       }
00566       
00567       /* and do it */
00568       if (!solve_reduced(n, m, h_x, h_y, a, delta_x, delta_y, c_x, c_y, workspace, PREDICTOR))
00569       {
00570           status=INCONSISTENT;
00571           goto exit_optimizer;
00572       }
00573       
00574       for (i=0; i<n; i++) {
00575     /* backsubstitution */
00576     delta_s[i] = s[i] * (delta_x[i] - hat_tau[i]) / t[i];
00577     delta_z[i] = z[i] * (hat_nu[i] - delta_x[i]) / g[i];
00578     
00579     delta_g[i] = g[i] * (gamma_z[i] - delta_z[i]) / z[i];
00580     delta_t[i] = t[i] * (gamma_s[i] - delta_s[i]) / s[i];
00581     
00582     /* central path (corrector) */
00583     gamma_z[i] = mu / g[i] - z[i] - delta_z[i] * delta_g[i] / g[i];
00584     gamma_s[i] = mu / t[i] - s[i] - delta_s[i] * delta_t[i] / t[i];
00585     
00586     /* (some more intermediate variables) the hat variables */
00587     hat_nu[i] = nu[i] + g[i] * gamma_z[i] / z[i];
00588     hat_tau[i] = tau[i] - t[i] * gamma_s[i] / s[i];
00589     
00590     /* initialization before the cholesky */
00591     c_x[i] = sigma[i] - z[i] * hat_nu[i] / g[i] - s[i] * hat_tau[i] / t[i];
00592       }
00593       
00594       for (i=0; i<m; i++) { /* comput c_y and rho */
00595     c_y[i] = rho[i];
00596     for (j=i; j<m; j++)
00597       h_y[m*i + j] = 0;
00598       }
00599       
00600       /* and do it */
00601       solve_reduced(n, m, h_x, h_y, a, delta_x, delta_y, c_x, c_y, workspace,
00602             CORRECTOR);
00603       
00604       for (i=0; i<n; i++) {
00605     /* backsubstitution */
00606     delta_s[i] = s[i] * (delta_x[i] - hat_tau[i]) / t[i];
00607     delta_z[i] = z[i] * (hat_nu[i] - delta_x[i]) / g[i];
00608     
00609     delta_g[i] = g[i] * (gamma_z[i] - delta_z[i]) / z[i];
00610     delta_t[i] = t[i] * (gamma_s[i] - delta_s[i]) / s[i];
00611       }
00612       
00613       alfa = -1;
00614       for (i=0; i<n; i++) {
00615     alfa = CMath::min(alfa, delta_g[i]/g[i]);
00616     alfa = CMath::min(alfa, delta_t[i]/t[i]);
00617     alfa = CMath::min(alfa, delta_s[i]/s[i]);
00618     alfa = CMath::min(alfa, delta_z[i]/z[i]);
00619       }
00620       alfa = (margin - 1) / alfa;
00621       
00622       /* compute mu */
00623       for (i=0, mu=0; i<n; i++)
00624     mu += z[i] * g[i] + s[i] * t[i];
00625       mu = mu / (2*n);
00626       mu = mu * CMath::sq((alfa - 1) / (alfa + 10));
00627       
00628       for (i=0; i<n; i++) {
00629     x[i] += alfa * delta_x[i];
00630     g[i] += alfa * delta_g[i];
00631     t[i] += alfa * delta_t[i];
00632     z[i] += alfa * delta_z[i];
00633     s[i] += alfa * delta_s[i];
00634       }
00635       
00636       for (i=0; i<m; i++) 
00637     y[i] += alfa * delta_y[i];
00638     }
00639   }
00640 
00641 exit_optimizer: 
00642   if ((status == 1) && (verb >= STATUS)) {
00643       SG_SDEBUG("----------------------------------------------------------------------------------\n");
00644       SG_SDEBUG("optimization converged\n");
00645   }
00646   
00647   /* free memory */
00648   free(workspace);
00649   free(diag_h_x);
00650   free(h_y);
00651   free(c_x);
00652   free(c_y);
00653   free(h_dot_x);
00654   
00655   free(rho);
00656   free(nu);
00657   free(tau);
00658   free(sigma);
00659   free(gamma_z);
00660   free(gamma_s);
00661   
00662   free(hat_nu);
00663   free(hat_tau);
00664     
00665   free(delta_x);
00666   free(delta_y);
00667   free(delta_s);
00668   free(delta_z);
00669   free(delta_g);
00670   free(delta_t);
00671     
00672   free(d);
00673 
00674   /* and return to sender */
00675   return status;
00676 }

SHOGUN Machine Learning Toolbox - Documentation