From 1ade55b50efc9f64b2fd44c42a324cd4b05b24b6 Mon Sep 17 00:00:00 2001
From: Julien Malik <julien.malik@c-s.fr>
Date: Thu, 14 Oct 2010 13:12:14 +0200
Subject: [PATCH] ENH: update libsvm to version 3.0
---
Code/Learning/otbSVMModel.txx | 8 +-
...leKernelFunctorImageClassificationTest.cxx | 1 -
...leKernelFunctorImageModelEstimatorTest.cxx | 1 -
Utilities/otbsvm/CMakeLists.txt | 18 +-
Utilities/otbsvm/svm.cpp | 3790 +++++++++++++++++
Utilities/otbsvm/svm.cxx | 3788 ----------------
Utilities/otbsvm/svm.h | 425 +-
7 files changed, 4025 insertions(+), 4006 deletions(-)
create mode 100644 Utilities/otbsvm/svm.cpp
delete mode 100644 Utilities/otbsvm/svm.cxx
diff --git a/Code/Learning/otbSVMModel.txx b/Code/Learning/otbSVMModel.txx
index 6c70c282f7..18dc0e4c23 100644
--- a/Code/Learning/otbSVMModel.txx
+++ b/Code/Learning/otbSVMModel.txx
@@ -111,12 +111,10 @@ template <class TValue, class TLabel>
void
SVMModel<TValue, TLabel>::DeleteModel()
{
- if (m_Model)
- {
- svm_destroy_model(m_Model);
- m_Model = NULL;
- }
+ svm_free_and_destroy_model(&m_Model);
+ m_Model = NULL;
}
+
template <class TValue, class TLabel>
void
SVMModel<TValue, TLabel>::DeleteProblem()
diff --git a/Testing/Code/Learning/otbSVMInverseCosSpectralAngleKernelFunctorImageClassificationTest.cxx b/Testing/Code/Learning/otbSVMInverseCosSpectralAngleKernelFunctorImageClassificationTest.cxx
index 0f06206df1..0de7effa0a 100644
--- a/Testing/Code/Learning/otbSVMInverseCosSpectralAngleKernelFunctorImageClassificationTest.cxx
+++ b/Testing/Code/Learning/otbSVMInverseCosSpectralAngleKernelFunctorImageClassificationTest.cxx
@@ -64,7 +64,6 @@ int otbSVMInverseCosSpectralAngleKernelFunctorImageClassificationTest(int argc,
otb::InverseCosSAMKernelFunctor myKernel;
myKernel.SetValue("Coef", 1.0);
- myKernel.Update();
svmEstimator->SetKernelFunctor(&myKernel);
svmEstimator->SetKernelType(GENERIC);
diff --git a/Testing/Code/Learning/otbSVMInverseCosSpectralAngleKernelFunctorImageModelEstimatorTest.cxx b/Testing/Code/Learning/otbSVMInverseCosSpectralAngleKernelFunctorImageModelEstimatorTest.cxx
index c63c5c4d99..2b353ecc18 100644
--- a/Testing/Code/Learning/otbSVMInverseCosSpectralAngleKernelFunctorImageModelEstimatorTest.cxx
+++ b/Testing/Code/Learning/otbSVMInverseCosSpectralAngleKernelFunctorImageModelEstimatorTest.cxx
@@ -63,7 +63,6 @@ int otbSVMInverseCosSpectralAngleKernelFunctorImageModelEstimatorTest(int argc,
otb::InverseCosSAMKernelFunctor myKernel;
myKernel.SetValue("Coef", 1.0);
- myKernel.Update();
svmEstimator->SetKernelFunctor(&myKernel);
svmEstimator->SetKernelType(GENERIC);
diff --git a/Utilities/otbsvm/CMakeLists.txt b/Utilities/otbsvm/CMakeLists.txt
index d517b80355..62591f3e9c 100644
--- a/Utilities/otbsvm/CMakeLists.txt
+++ b/Utilities/otbsvm/CMakeLists.txt
@@ -8,25 +8,9 @@ PROJECT(OTBSVM)
# source files for otbsvm
SET(OTBSVM_SRCS
-svm.cxx
+svm.cpp
)
-# Adjust the compiler flags to avoid problems with ossim code.
-IF(CMAKE_COMPILER_IS_GNUCXX)
- # Hide warnings in ossim code.
- SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -w")
-
- # CMake versions before 2.4 did not put the source-file-specific flags
- # after the directory-level flags. Remove optimization flags from the
- # CMAKE_CXX_FLAGS* variables so they can be overridden.
- IF("${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}" LESS 2.4)
- FOREACH(c "" "_DEBUG" "_RELEASE" "_MINSIZEREL" "_RELWITHDEBINFO")
- STRING(REGEX REPLACE
- "-O[^ ]*( |$)" "" CMAKE_CXX_FLAGS${c} "${CMAKE_CXX_FLAGS${c}}")
- ENDFOREACH(c)
- ENDIF("${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}" LESS 2.4)
-ENDIF(CMAKE_COMPILER_IS_GNUCXX)
-
ADD_LIBRARY(otbsvm ${OTBSVM_SRCS})
TARGET_LINK_LIBRARIES (otbsvm ITKCommon)
IF(OTB_LIBRARY_PROPERTIES)
diff --git a/Utilities/otbsvm/svm.cpp b/Utilities/otbsvm/svm.cpp
new file mode 100644
index 0000000000..d566a9c964
--- /dev/null
+++ b/Utilities/otbsvm/svm.cpp
@@ -0,0 +1,3790 @@
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <ctype.h>
+#include <float.h>
+#include <string.h>
+#include <stdarg.h>
+#include "svm.h"
+int libsvm_version = LIBSVM_VERSION;
+typedef float Qfloat;
+typedef signed char schar;
+#ifndef min
+template <class T> static inline T min(T x,T y) { return (x<y)?x:y; }
+#endif
+#ifndef max
+template <class T> static inline T max(T x,T y) { return (x>y)?x:y; }
+#endif
+template <class T> static inline void swap(T& x, T& y) { T t=x; x=y; y=t; }
+template <class S, class T> static inline void clone(T*& dst, S* src, int n)
+{
+ dst = new T[n];
+ memcpy((void *)dst,(void *)src,sizeof(T)*n);
+}
+static inline double powi(double base, int times)
+{
+ double tmp = base, ret = 1.0;
+
+ for(int t=times; t>0; t/=2)
+ {
+ if(t%2==1) ret*=tmp;
+ tmp = tmp * tmp;
+ }
+ return ret;
+}
+#define INF HUGE_VAL
+#define TAU 1e-12
+#define Malloc(type,n) (type *)malloc((n)*sizeof(type))
+
+static void print_string_stdout(const char *s)
+{
+ fputs(s,stdout);
+ fflush(stdout);
+}
+static void (*svm_print_string) (const char *) = &print_string_stdout;
+#if 1
+static void info(const char *fmt,...)
+{
+ char buf[BUFSIZ];
+ va_list ap;
+ va_start(ap,fmt);
+ vsprintf(buf,fmt,ap);
+ va_end(ap);
+ (*svm_print_string)(buf);
+}
+#else
+static void info(const char *fmt,...) {}
+#endif
+
+//
+// Kernel Cache
+//
+// l is the number of total data items
+// size is the cache size limit in bytes
+//
+class Cache
+{
+public:
+ Cache(int l,long int size);
+ ~Cache();
+
+ // request data [0,len)
+ // return some position p where [p,len) need to be filled
+ // (p >= len if nothing needs to be filled)
+ int get_data(const int index, Qfloat **data, int len);
+ void swap_index(int i, int j);
+private:
+ int l;
+ long int size;
+ struct head_t
+ {
+ head_t *prev, *next; // a circular list
+ Qfloat *data;
+ int len; // data[0,len) is cached in this entry
+ };
+
+ head_t *head;
+ head_t lru_head;
+ void lru_delete(head_t *h);
+ void lru_insert(head_t *h);
+};
+
+Cache::Cache(int l_,long int size_):l(l_),size(size_)
+{
+ head = (head_t *)calloc(l,sizeof(head_t)); // initialized to 0
+ size /= sizeof(Qfloat);
+ size -= l * sizeof(head_t) / sizeof(Qfloat);
+ size = max(size, 2 * (long int) l); // cache must be large enough for two columns
+ lru_head.next = lru_head.prev = &lru_head;
+}
+
+Cache::~Cache()
+{
+ for(head_t *h = lru_head.next; h != &lru_head; h=h->next)
+ free(h->data);
+ free(head);
+}
+
+void Cache::lru_delete(head_t *h)
+{
+ // delete from current location
+ h->prev->next = h->next;
+ h->next->prev = h->prev;
+}
+
+void Cache::lru_insert(head_t *h)
+{
+ // insert to last position
+ h->next = &lru_head;
+ h->prev = lru_head.prev;
+ h->prev->next = h;
+ h->next->prev = h;
+}
+
+int Cache::get_data(const int index, Qfloat **data, int len)
+{
+ head_t *h = &head[index];
+ if(h->len) lru_delete(h);
+ int more = len - h->len;
+
+ if(more > 0)
+ {
+ // free old space
+ while(size < more)
+ {
+ head_t *old = lru_head.next;
+ lru_delete(old);
+ free(old->data);
+ size += old->len;
+ old->data = 0;
+ old->len = 0;
+ }
+
+ // allocate new space
+ h->data = (Qfloat *)realloc(h->data,sizeof(Qfloat)*len);
+ size -= more;
+ swap(h->len,len);
+ }
+
+ lru_insert(h);
+ *data = h->data;
+ return len;
+}
+
+void Cache::swap_index(int i, int j)
+{
+ if(i==j) return;
+
+ if(head[i].len) lru_delete(&head[i]);
+ if(head[j].len) lru_delete(&head[j]);
+ swap(head[i].data,head[j].data);
+ swap(head[i].len,head[j].len);
+ if(head[i].len) lru_insert(&head[i]);
+ if(head[j].len) lru_insert(&head[j]);
+
+ if(i>j) swap(i,j);
+ for(head_t *h = lru_head.next; h!=&lru_head; h=h->next)
+ {
+ if(h->len > i)
+ {
+ if(h->len > j)
+ swap(h->data[i],h->data[j]);
+ else
+ {
+ // give up
+ lru_delete(h);
+ free(h->data);
+ size += h->len;
+ h->data = 0;
+ h->len = 0;
+ }
+ }
+ }
+}
+
+//
+// Kernel evaluation
+//
+// the static method k_function is for doing single kernel evaluation
+// the constructor of Kernel prepares to calculate the l*l kernel matrix
+// the member function get_Q is for getting one column from the Q Matrix
+//
+class QMatrix {
+public:
+ virtual Qfloat *get_Q(int column, int len) const = 0;
+ virtual double *get_QD() const = 0;
+ virtual void swap_index(int i, int j) const = 0;
+ virtual ~QMatrix() {}
+};
+
+class Kernel: public QMatrix {
+public:
+ Kernel(int l, svm_node * const * x, const svm_parameter& param);
+ virtual ~Kernel();
+
+ static double k_function(const svm_node *x, const svm_node *y,
+ const svm_parameter& param);
+ virtual Qfloat *get_Q(int column, int len) const = 0;
+ virtual double *get_QD() const = 0;
+ virtual void swap_index(int i, int j) const // no so const...
+ {
+ swap(x[i],x[j]);
+ if(x_square) swap(x_square[i],x_square[j]);
+ }
+protected:
+
+ /*** Begin OTB modification ***/
+ //double (Kernel::*kernel_function)(int i, int j) const;
+ double (Kernel::*kernel_function)(int i, int j, const svm_parameter& param) const;
+ const svm_parameter& m_param;
+ /*** End OTB modification ***/
+
+private:
+ const svm_node **x;
+ double *x_square;
+
+ // svm_parameter
+ const int kernel_type;
+ const int degree;
+ const double gamma;
+ const double coef0;
+
+ static double dot(const svm_node *px, const svm_node *py);
+
+ /*** OTB modification : add svm_parameter to the list of parameters ***/
+ double kernel_linear(int i, int j, const svm_parameter& param) const
+ {
+ return dot(x[i],x[j]);
+ }
+ double kernel_poly(int i, int j, const svm_parameter& param) const
+ {
+ return powi(gamma*dot(x[i],x[j])+coef0,degree);
+ }
+ double kernel_rbf(int i, int j, const svm_parameter& param) const
+ {
+ return exp(-gamma*(x_square[i]+x_square[j]-2*dot(x[i],x[j])));
+ }
+ double kernel_sigmoid(int i, int j, const svm_parameter& param) const
+ {
+ return tanh(gamma*dot(x[i],x[j])+coef0);
+ }
+ double kernel_precomputed(int i, int j, const svm_parameter& param) const
+ {
+ return x[i][(int)(x[j][0].value)].value;
+ }
+
+ /*** Begin OTB modification ***/
+ double kernel_generic(int i, int j, const svm_parameter& param) const
+ {
+ if( param.kernel_generic == NULL )
+ {
+ itkGenericExceptionMacro( << "Generic Kernel is not initialized !");
+ }
+ return (*param.kernel_generic)(x[i],x[j],param);
+ }
+
+ double kernel_composed(int i, int j, const svm_parameter& param) const
+ {
+ if( param.kernel_composed == NULL )
+ {
+ itkGenericExceptionMacro( << "Generic Kernel is not initialized !");
+ }
+ return (*param.kernel_composed)(x[i],x[j],param);
+ }
+ /*** End OTB modification ***/
+};
+
+Kernel::Kernel(int l, svm_node * const * x_, const svm_parameter& param)
+:kernel_type(param.kernel_type), degree(param.degree),
+ gamma(param.gamma), coef0(param.coef0)
+/*** Begin OTB modification ***/
+, m_param(param)
+/*** End OTB modification ***/
+{
+ switch(kernel_type)
+ {
+ case LINEAR:
+ kernel_function = &Kernel::kernel_linear;
+ break;
+ case POLY:
+ kernel_function = &Kernel::kernel_poly;
+ break;
+ case RBF:
+ kernel_function = &Kernel::kernel_rbf;
+ break;
+ case SIGMOID:
+ kernel_function = &Kernel::kernel_sigmoid;
+ break;
+ case PRECOMPUTED:
+ kernel_function = &Kernel::kernel_precomputed;
+ break;
+ /*** Begin OTB modification ***/
+ case GENERIC:
+ kernel_function = &Kernel::kernel_generic;
+ break;
+ case COMPOSED:
+ kernel_function = &Kernel::kernel_composed;
+ break;
+ /*** End OTB modification ***/
+ }
+
+ clone(x,x_,l);
+
+ if(kernel_type == RBF)
+ {
+ x_square = new double[l];
+ for(int i=0;i<l;i++)
+ x_square[i] = dot(x[i],x[i]);
+ }
+ else
+ x_square = 0;
+}
+
+Kernel::~Kernel()
+{
+ delete[] x;
+ delete[] x_square;
+}
+
+double Kernel::dot(const svm_node *px, const svm_node *py)
+{
+ double sum = 0;
+ while(px->index != -1 && py->index != -1)
+ {
+ if(px->index == py->index)
+ {
+ sum += px->value * py->value;
+ ++px;
+ ++py;
+ }
+ else
+ {
+ if(px->index > py->index)
+ ++py;
+ else
+ ++px;
+ }
+ }
+ return sum;
+}
+
+double Kernel::k_function(const svm_node *x, const svm_node *y,
+ const svm_parameter& param)
+{
+ switch(param.kernel_type)
+ {
+ case LINEAR:
+ return dot(x,y);
+ case POLY:
+ return powi(param.gamma*dot(x,y)+param.coef0,param.degree);
+ case RBF:
+ {
+ double sum = 0;
+ while(x->index != -1 && y->index !=-1)
+ {
+ if(x->index == y->index)
+ {
+ double d = x->value - y->value;
+ sum += d*d;
+ ++x;
+ ++y;
+ }
+ else
+ {
+ if(x->index > y->index)
+ {
+ sum += y->value * y->value;
+ ++y;
+ }
+ else
+ {
+ sum += x->value * x->value;
+ ++x;
+ }
+ }
+ }
+
+ while(x->index != -1)
+ {
+ sum += x->value * x->value;
+ ++x;
+ }
+
+ while(y->index != -1)
+ {
+ sum += y->value * y->value;
+ ++y;
+ }
+
+ return exp(-param.gamma*sum);
+ }
+ case SIGMOID:
+ return tanh(param.gamma*dot(x,y)+param.coef0);
+ case PRECOMPUTED: //x: test (validation), y: SV
+ return x[(int)(y->value)].value;
+
+ /*** Begin OTB modification ***/
+ case GENERIC:
+ return ((*param.kernel_generic)(x,y,param));
+ case COMPOSED:
+ return ((*param.kernel_composed)(x,y,param));
+ /*** End OTB modification ***/
+
+ default:
+ return 0; // Unreachable
+ }
+}
+
+// An SMO algorithm in Fan et al., JMLR 6(2005), p. 1889--1918
+// Solves:
+//
+// min 0.5(\alpha^T Q \alpha) + p^T \alpha
+//
+// y^T \alpha = \delta
+// y_i = +1 or -1
+// 0 <= alpha_i <= Cp for y_i = 1
+// 0 <= alpha_i <= Cn for y_i = -1
+//
+// Given:
+//
+// Q, p, y, Cp, Cn, and an initial feasible point \alpha
+// l is the size of vectors and matrices
+// eps is the stopping tolerance
+//
+// solution will be put in \alpha, objective value will be put in obj
+//
+class Solver {
+public:
+ Solver() {};
+ virtual ~Solver() {};
+
+ struct SolutionInfo {
+ double obj;
+ double rho;
+ double upper_bound_p;
+ double upper_bound_n;
+ double r; // for Solver_NU
+ };
+
+ void Solve(int l, const QMatrix& Q, const double *p_, const schar *y_,
+ double *alpha_, double Cp, double Cn, double eps,
+ SolutionInfo* si, int shrinking);
+protected:
+ int active_size;
+ schar *y;
+ double *G; // gradient of objective function
+ enum { LOWER_BOUND, UPPER_BOUND, FREE };
+ char *alpha_status; // LOWER_BOUND, UPPER_BOUND, FREE
+ double *alpha;
+ const QMatrix *Q;
+ const double *QD;
+ double eps;
+ double Cp,Cn;
+ double *p;
+ int *active_set;
+ double *G_bar; // gradient, if we treat free variables as 0
+ int l;
+ bool unshrink; // XXX
+
+ double get_C(int i)
+ {
+ return (y[i] > 0)? Cp : Cn;
+ }
+ void update_alpha_status(int i)
+ {
+ if(alpha[i] >= get_C(i))
+ alpha_status[i] = UPPER_BOUND;
+ else if(alpha[i] <= 0)
+ alpha_status[i] = LOWER_BOUND;
+ else alpha_status[i] = FREE;
+ }
+ bool is_upper_bound(int i) { return alpha_status[i] == UPPER_BOUND; }
+ bool is_lower_bound(int i) { return alpha_status[i] == LOWER_BOUND; }
+ bool is_free(int i) { return alpha_status[i] == FREE; }
+ void swap_index(int i, int j);
+ void reconstruct_gradient();
+ virtual int select_working_set(int &i, int &j);
+ virtual double calculate_rho();
+ virtual void do_shrinking();
+private:
+ bool be_shrunk(int i, double Gmax1, double Gmax2);
+};
+
+void Solver::swap_index(int i, int j)
+{
+ Q->swap_index(i,j);
+ swap(y[i],y[j]);
+ swap(G[i],G[j]);
+ swap(alpha_status[i],alpha_status[j]);
+ swap(alpha[i],alpha[j]);
+ swap(p[i],p[j]);
+ swap(active_set[i],active_set[j]);
+ swap(G_bar[i],G_bar[j]);
+}
+
+void Solver::reconstruct_gradient()
+{
+ // reconstruct inactive elements of G from G_bar and free variables
+
+ if(active_size == l) return;
+
+ int i,j;
+ int nr_free = 0;
+
+ for(j=active_size;j<l;j++)
+ G[j] = G_bar[j] + p[j];
+
+ for(j=0;j<active_size;j++)
+ if(is_free(j))
+ nr_free++;
+
+ if(2*nr_free < active_size)
+ info("\nWarning: using -h 0 may be faster\n");
+
+ if (nr_free*l > 2*active_size*(l-active_size))
+ {
+ for(i=active_size;i<l;i++)
+ {
+ const Qfloat *Q_i = Q->get_Q(i,active_size);
+ for(j=0;j<active_size;j++)
+ if(is_free(j))
+ G[i] += alpha[j] * Q_i[j];
+ }
+ }
+ else
+ {
+ for(i=0;i<active_size;i++)
+ if(is_free(i))
+ {
+ const Qfloat *Q_i = Q->get_Q(i,l);
+ double alpha_i = alpha[i];
+ for(j=active_size;j<l;j++)
+ G[j] += alpha_i * Q_i[j];
+ }
+ }
+}
+
+void Solver::Solve(int l, const QMatrix& Q, const double *p_, const schar *y_,
+ double *alpha_, double Cp, double Cn, double eps,
+ SolutionInfo* si, int shrinking)
+{
+ this->l = l;
+ this->Q = &Q;
+ QD=Q.get_QD();
+ clone(p, p_,l);
+ clone(y, y_,l);
+ clone(alpha,alpha_,l);
+ this->Cp = Cp;
+ this->Cn = Cn;
+ this->eps = eps;
+ unshrink = false;
+
+ // initialize alpha_status
+ {
+ alpha_status = new char[l];
+ for(int i=0;i<l;i++)
+ update_alpha_status(i);
+ }
+
+ // initialize active set (for shrinking)
+ {
+ active_set = new int[l];
+ for(int i=0;i<l;i++)
+ active_set[i] = i;
+ active_size = l;
+ }
+
+ // initialize gradient
+ {
+ G = new double[l];
+ G_bar = new double[l];
+ int i;
+ for(i=0;i<l;i++)
+ {
+ G[i] = p[i];
+ G_bar[i] = 0;
+ }
+ for(i=0;i<l;i++)
+ if(!is_lower_bound(i))
+ {
+ const Qfloat *Q_i = Q.get_Q(i,l);
+ double alpha_i = alpha[i];
+ int j;
+ for(j=0;j<l;j++)
+ G[j] += alpha_i*Q_i[j];
+ if(is_upper_bound(i))
+ for(j=0;j<l;j++)
+ G_bar[j] += get_C(i) * Q_i[j];
+ }
+ }
+
+ // optimization step
+
+ int iter = 0;
+ int counter = min(l,1000)+1;
+
+ while(1)
+ {
+ // show progress and do shrinking
+
+ if(--counter == 0)
+ {
+ counter = min(l,1000);
+ if(shrinking) do_shrinking();
+ info(".");
+ }
+
+ int i,j;
+ if(select_working_set(i,j)!=0)
+ {
+ // reconstruct the whole gradient
+ reconstruct_gradient();
+ // reset active set size and check
+ active_size = l;
+ info("*");
+ if(select_working_set(i,j)!=0)
+ break;
+ else
+ counter = 1; // do shrinking next iteration
+ }
+
+ ++iter;
+
+ // update alpha[i] and alpha[j], handle bounds carefully
+
+ const Qfloat *Q_i = Q.get_Q(i,active_size);
+ const Qfloat *Q_j = Q.get_Q(j,active_size);
+
+ double C_i = get_C(i);
+ double C_j = get_C(j);
+
+ double old_alpha_i = alpha[i];
+ double old_alpha_j = alpha[j];
+
+ if(y[i]!=y[j])
+ {
+ double quad_coef = QD[i]+QD[j]+2*Q_i[j];
+ if (quad_coef <= 0)
+ quad_coef = TAU;
+ double delta = (-G[i]-G[j])/quad_coef;
+ double diff = alpha[i] - alpha[j];
+ alpha[i] += delta;
+ alpha[j] += delta;
+
+ if(diff > 0)
+ {
+ if(alpha[j] < 0)
+ {
+ alpha[j] = 0;
+ alpha[i] = diff;
+ }
+ }
+ else
+ {
+ if(alpha[i] < 0)
+ {
+ alpha[i] = 0;
+ alpha[j] = -diff;
+ }
+ }
+ if(diff > C_i - C_j)
+ {
+ if(alpha[i] > C_i)
+ {
+ alpha[i] = C_i;
+ alpha[j] = C_i - diff;
+ }
+ }
+ else
+ {
+ if(alpha[j] > C_j)
+ {
+ alpha[j] = C_j;
+ alpha[i] = C_j + diff;
+ }
+ }
+ }
+ else
+ {
+ double quad_coef = QD[i]+QD[j]-2*Q_i[j];
+ if (quad_coef <= 0)
+ quad_coef = TAU;
+ double delta = (G[i]-G[j])/quad_coef;
+ double sum = alpha[i] + alpha[j];
+ alpha[i] -= delta;
+ alpha[j] += delta;
+
+ if(sum > C_i)
+ {
+ if(alpha[i] > C_i)
+ {
+ alpha[i] = C_i;
+ alpha[j] = sum - C_i;
+ }
+ }
+ else
+ {
+ if(alpha[j] < 0)
+ {
+ alpha[j] = 0;
+ alpha[i] = sum;
+ }
+ }
+ if(sum > C_j)
+ {
+ if(alpha[j] > C_j)
+ {
+ alpha[j] = C_j;
+ alpha[i] = sum - C_j;
+ }
+ }
+ else
+ {
+ if(alpha[i] < 0)
+ {
+ alpha[i] = 0;
+ alpha[j] = sum;
+ }
+ }
+ }
+
+ // update G
+
+ double delta_alpha_i = alpha[i] - old_alpha_i;
+ double delta_alpha_j = alpha[j] - old_alpha_j;
+
+ for(int k=0;k<active_size;k++)
+ {
+ G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j;
+ }
+
+ // update alpha_status and G_bar
+
+ {
+ bool ui = is_upper_bound(i);
+ bool uj = is_upper_bound(j);
+ update_alpha_status(i);
+ update_alpha_status(j);
+ int k;
+ if(ui != is_upper_bound(i))
+ {
+ Q_i = Q.get_Q(i,l);
+ if(ui)
+ for(k=0;k<l;k++)
+ G_bar[k] -= C_i * Q_i[k];
+ else
+ for(k=0;k<l;k++)
+ G_bar[k] += C_i * Q_i[k];
+ }
+
+ if(uj != is_upper_bound(j))
+ {
+ Q_j = Q.get_Q(j,l);
+ if(uj)
+ for(k=0;k<l;k++)
+ G_bar[k] -= C_j * Q_j[k];
+ else
+ for(k=0;k<l;k++)
+ G_bar[k] += C_j * Q_j[k];
+ }
+ }
+ }
+
+ // calculate rho
+
+ si->rho = calculate_rho();
+
+ // calculate objective value
+ {
+ double v = 0;
+ int i;
+ for(i=0;i<l;i++)
+ v += alpha[i] * (G[i] + p[i]);
+
+ si->obj = v/2;
+ }
+
+ // put back the solution
+ {
+ for(int i=0;i<l;i++)
+ alpha_[active_set[i]] = alpha[i];
+ }
+
+ // juggle everything back
+ /*{
+ for(int i=0;i<l;i++)
+ while(active_set[i] != i)
+ swap_index(i,active_set[i]);
+ // or Q.swap_index(i,active_set[i]);
+ }*/
+
+ si->upper_bound_p = Cp;
+ si->upper_bound_n = Cn;
+
+ info("\noptimization finished, #iter = %d\n",iter);
+
+ delete[] p;
+ delete[] y;
+ delete[] alpha;
+ delete[] alpha_status;
+ delete[] active_set;
+ delete[] G;
+ delete[] G_bar;
+}
+
+// return 1 if already optimal, return 0 otherwise
+int Solver::select_working_set(int &out_i, int &out_j)
+{
+ // return i,j such that
+ // i: maximizes -y_i * grad(f)_i, i in I_up(\alpha)
+ // j: minimizes the decrease of obj value
+ // (if quadratic coefficeint <= 0, replace it with tau)
+ // -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha)
+
+ double Gmax = -INF;
+ double Gmax2 = -INF;
+ int Gmax_idx = -1;
+ int Gmin_idx = -1;
+ double obj_diff_min = INF;
+
+ for(int t=0;t<active_size;t++)
+ if(y[t]==+1)
+ {
+ if(!is_upper_bound(t))
+ if(-G[t] >= Gmax)
+ {
+ Gmax = -G[t];
+ Gmax_idx = t;
+ }
+ }
+ else
+ {
+ if(!is_lower_bound(t))
+ if(G[t] >= Gmax)
+ {
+ Gmax = G[t];
+ Gmax_idx = t;
+ }
+ }
+
+ int i = Gmax_idx;
+ const Qfloat *Q_i = NULL;
+ if(i != -1) // NULL Q_i not accessed: Gmax=-INF if i=-1
+ Q_i = Q->get_Q(i,active_size);
+
+ for(int j=0;j<active_size;j++)
+ {
+ if(y[j]==+1)
+ {
+ if (!is_lower_bound(j))
+ {
+ double grad_diff=Gmax+G[j];
+ if (G[j] >= Gmax2)
+ Gmax2 = G[j];
+ if (grad_diff > 0)
+ {
+ double obj_diff;
+ double quad_coef = QD[i]+QD[j]-2.0*y[i]*Q_i[j];
+ if (quad_coef > 0)
+ obj_diff = -(grad_diff*grad_diff)/quad_coef;
+ else
+ obj_diff = -(grad_diff*grad_diff)/TAU;
+
+ if (obj_diff <= obj_diff_min)
+ {
+ Gmin_idx=j;
+ obj_diff_min = obj_diff;
+ }
+ }
+ }
+ }
+ else
+ {
+ if (!is_upper_bound(j))
+ {
+ double grad_diff= Gmax-G[j];
+ if (-G[j] >= Gmax2)
+ Gmax2 = -G[j];
+ if (grad_diff > 0)
+ {
+ double obj_diff;
+ double quad_coef = QD[i]+QD[j]+2.0*y[i]*Q_i[j];
+ if (quad_coef > 0)
+ obj_diff = -(grad_diff*grad_diff)/quad_coef;
+ else
+ obj_diff = -(grad_diff*grad_diff)/TAU;
+
+ if (obj_diff <= obj_diff_min)
+ {
+ Gmin_idx=j;
+ obj_diff_min = obj_diff;
+ }
+ }
+ }
+ }
+ }
+
+ if(Gmax+Gmax2 < eps)
+ return 1;
+
+ out_i = Gmax_idx;
+ out_j = Gmin_idx;
+ return 0;
+}
+
+bool Solver::be_shrunk(int i, double Gmax1, double Gmax2)
+{
+ if(is_upper_bound(i))
+ {
+ if(y[i]==+1)
+ return(-G[i] > Gmax1);
+ else
+ return(-G[i] > Gmax2);
+ }
+ else if(is_lower_bound(i))
+ {
+ if(y[i]==+1)
+ return(G[i] > Gmax2);
+ else
+ return(G[i] > Gmax1);
+ }
+ else
+ return(false);
+}
+
+void Solver::do_shrinking()
+{
+ int i;
+ double Gmax1 = -INF; // max { -y_i * grad(f)_i | i in I_up(\alpha) }
+ double Gmax2 = -INF; // max { y_i * grad(f)_i | i in I_low(\alpha) }
+
+ // find maximal violating pair first
+ for(i=0;i<active_size;i++)
+ {
+ if(y[i]==+1)
+ {
+ if(!is_upper_bound(i))
+ {
+ if(-G[i] >= Gmax1)
+ Gmax1 = -G[i];
+ }
+ if(!is_lower_bound(i))
+ {
+ if(G[i] >= Gmax2)
+ Gmax2 = G[i];
+ }
+ }
+ else
+ {
+ if(!is_upper_bound(i))
+ {
+ if(-G[i] >= Gmax2)
+ Gmax2 = -G[i];
+ }
+ if(!is_lower_bound(i))
+ {
+ if(G[i] >= Gmax1)
+ Gmax1 = G[i];
+ }
+ }
+ }
+
+ if(unshrink == false && Gmax1 + Gmax2 <= eps*10)
+ {
+ unshrink = true;
+ reconstruct_gradient();
+ active_size = l;
+ info("*");
+ }
+
+ for(i=0;i<active_size;i++)
+ if (be_shrunk(i, Gmax1, Gmax2))
+ {
+ active_size--;
+ while (active_size > i)
+ {
+ if (!be_shrunk(active_size, Gmax1, Gmax2))
+ {
+ swap_index(i,active_size);
+ break;
+ }
+ active_size--;
+ }
+ }
+}
+
+double Solver::calculate_rho()
+{
+ double r;
+ int nr_free = 0;
+ double ub = INF, lb = -INF, sum_free = 0;
+ for(int i=0;i<active_size;i++)
+ {
+ double yG = y[i]*G[i];
+
+ if(is_upper_bound(i))
+ {
+ if(y[i]==-1)
+ ub = min(ub,yG);
+ else
+ lb = max(lb,yG);
+ }
+ else if(is_lower_bound(i))
+ {
+ if(y[i]==+1)
+ ub = min(ub,yG);
+ else
+ lb = max(lb,yG);
+ }
+ else
+ {
+ ++nr_free;
+ sum_free += yG;
+ }
+ }
+
+ if(nr_free>0)
+ r = sum_free/nr_free;
+ else
+ r = (ub+lb)/2;
+
+ return r;
+}
+
+//
+// Solver for nu-svm classification and regression
+//
+// additional constraint: e^T \alpha = constant
+//
+class Solver_NU : public Solver
+{
+public:
+ Solver_NU() {}
+ void Solve(int l, const QMatrix& Q, const double *p, const schar *y,
+ double *alpha, double Cp, double Cn, double eps,
+ SolutionInfo* si, int shrinking)
+ {
+ this->si = si;
+ Solver::Solve(l,Q,p,y,alpha,Cp,Cn,eps,si,shrinking);
+ }
+private:
+ SolutionInfo *si;
+ int select_working_set(int &i, int &j);
+ double calculate_rho();
+ bool be_shrunk(int i, double Gmax1, double Gmax2, double Gmax3, double Gmax4);
+ void do_shrinking();
+};
+
+// return 1 if already optimal, return 0 otherwise
+int Solver_NU::select_working_set(int &out_i, int &out_j)
+{
+ // return i,j such that y_i = y_j and
+ // i: maximizes -y_i * grad(f)_i, i in I_up(\alpha)
+ // j: minimizes the decrease of obj value
+ // (if quadratic coefficeint <= 0, replace it with tau)
+ // -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha)
+
+ double Gmaxp = -INF;
+ double Gmaxp2 = -INF;
+ int Gmaxp_idx = -1;
+
+ double Gmaxn = -INF;
+ double Gmaxn2 = -INF;
+ int Gmaxn_idx = -1;
+
+ int Gmin_idx = -1;
+ double obj_diff_min = INF;
+
+ for(int t=0;t<active_size;t++)
+ if(y[t]==+1)
+ {
+ if(!is_upper_bound(t))
+ if(-G[t] >= Gmaxp)
+ {
+ Gmaxp = -G[t];
+ Gmaxp_idx = t;
+ }
+ }
+ else
+ {
+ if(!is_lower_bound(t))
+ if(G[t] >= Gmaxn)
+ {
+ Gmaxn = G[t];
+ Gmaxn_idx = t;
+ }
+ }
+
+ int ip = Gmaxp_idx;
+ int in = Gmaxn_idx;
+ const Qfloat *Q_ip = NULL;
+ const Qfloat *Q_in = NULL;
+ if(ip != -1) // NULL Q_ip not accessed: Gmaxp=-INF if ip=-1
+ Q_ip = Q->get_Q(ip,active_size);
+ if(in != -1)
+ Q_in = Q->get_Q(in,active_size);
+
+ for(int j=0;j<active_size;j++)
+ {
+ if(y[j]==+1)
+ {
+ if (!is_lower_bound(j))
+ {
+ double grad_diff=Gmaxp+G[j];
+ if (G[j] >= Gmaxp2)
+ Gmaxp2 = G[j];
+ if (grad_diff > 0)
+ {
+ double obj_diff;
+ double quad_coef = QD[ip]+QD[j]-2*Q_ip[j];
+ if (quad_coef > 0)
+ obj_diff = -(grad_diff*grad_diff)/quad_coef;
+ else
+ obj_diff = -(grad_diff*grad_diff)/TAU;
+
+ if (obj_diff <= obj_diff_min)
+ {
+ Gmin_idx=j;
+ obj_diff_min = obj_diff;
+ }
+ }
+ }
+ }
+ else
+ {
+ if (!is_upper_bound(j))
+ {
+ double grad_diff=Gmaxn-G[j];
+ if (-G[j] >= Gmaxn2)
+ Gmaxn2 = -G[j];
+ if (grad_diff > 0)
+ {
+ double obj_diff;
+ double quad_coef = QD[in]+QD[j]-2*Q_in[j];
+ if (quad_coef > 0)
+ obj_diff = -(grad_diff*grad_diff)/quad_coef;
+ else
+ obj_diff = -(grad_diff*grad_diff)/TAU;
+
+ if (obj_diff <= obj_diff_min)
+ {
+ Gmin_idx=j;
+ obj_diff_min = obj_diff;
+ }
+ }
+ }
+ }
+ }
+
+ if(max(Gmaxp+Gmaxp2,Gmaxn+Gmaxn2) < eps)
+ return 1;
+
+ if (y[Gmin_idx] == +1)
+ out_i = Gmaxp_idx;
+ else
+ out_i = Gmaxn_idx;
+ out_j = Gmin_idx;
+
+ return 0;
+}
+
+bool Solver_NU::be_shrunk(int i, double Gmax1, double Gmax2, double Gmax3, double Gmax4)
+{
+ if(is_upper_bound(i))
+ {
+ if(y[i]==+1)
+ return(-G[i] > Gmax1);
+ else
+ return(-G[i] > Gmax4);
+ }
+ else if(is_lower_bound(i))
+ {
+ if(y[i]==+1)
+ return(G[i] > Gmax2);
+ else
+ return(G[i] > Gmax3);
+ }
+ else
+ return(false);
+}
+
+void Solver_NU::do_shrinking()
+{
+ double Gmax1 = -INF; // max { -y_i * grad(f)_i | y_i = +1, i in I_up(\alpha) }
+ double Gmax2 = -INF; // max { y_i * grad(f)_i | y_i = +1, i in I_low(\alpha) }
+ double Gmax3 = -INF; // max { -y_i * grad(f)_i | y_i = -1, i in I_up(\alpha) }
+ double Gmax4 = -INF; // max { y_i * grad(f)_i | y_i = -1, i in I_low(\alpha) }
+
+ // find maximal violating pair first
+ int i;
+ for(i=0;i<active_size;i++)
+ {
+ if(!is_upper_bound(i))
+ {
+ if(y[i]==+1)
+ {
+ if(-G[i] > Gmax1) Gmax1 = -G[i];
+ }
+ else if(-G[i] > Gmax4) Gmax4 = -G[i];
+ }
+ if(!is_lower_bound(i))
+ {
+ if(y[i]==+1)
+ {
+ if(G[i] > Gmax2) Gmax2 = G[i];
+ }
+ else if(G[i] > Gmax3) Gmax3 = G[i];
+ }
+ }
+
+ if(unshrink == false && max(Gmax1+Gmax2,Gmax3+Gmax4) <= eps*10)
+ {
+ unshrink = true;
+ reconstruct_gradient();
+ active_size = l;
+ }
+
+ for(i=0;i<active_size;i++)
+ if (be_shrunk(i, Gmax1, Gmax2, Gmax3, Gmax4))
+ {
+ active_size--;
+ while (active_size > i)
+ {
+ if (!be_shrunk(active_size, Gmax1, Gmax2, Gmax3, Gmax4))
+ {
+ swap_index(i,active_size);
+ break;
+ }
+ active_size--;
+ }
+ }
+}
+
+double Solver_NU::calculate_rho()
+{
+ int nr_free1 = 0,nr_free2 = 0;
+ double ub1 = INF, ub2 = INF;
+ double lb1 = -INF, lb2 = -INF;
+ double sum_free1 = 0, sum_free2 = 0;
+
+ for(int i=0;i<active_size;i++)
+ {
+ if(y[i]==+1)
+ {
+ if(is_upper_bound(i))
+ lb1 = max(lb1,G[i]);
+ else if(is_lower_bound(i))
+ ub1 = min(ub1,G[i]);
+ else
+ {
+ ++nr_free1;
+ sum_free1 += G[i];
+ }
+ }
+ else
+ {
+ if(is_upper_bound(i))
+ lb2 = max(lb2,G[i]);
+ else if(is_lower_bound(i))
+ ub2 = min(ub2,G[i]);
+ else
+ {
+ ++nr_free2;
+ sum_free2 += G[i];
+ }
+ }
+ }
+
+ double r1,r2;
+ if(nr_free1 > 0)
+ r1 = sum_free1/nr_free1;
+ else
+ r1 = (ub1+lb1)/2;
+
+ if(nr_free2 > 0)
+ r2 = sum_free2/nr_free2;
+ else
+ r2 = (ub2+lb2)/2;
+
+ si->r = (r1+r2)/2;
+ return (r1-r2)/2;
+}
+
+//
+// Q matrices for various formulations
+//
+class SVC_Q: public Kernel
+{
+public:
+ SVC_Q(const svm_problem& prob, const svm_parameter& param, const schar *y_)
+ :Kernel(prob.l, prob.x, param)
+ {
+ clone(y,y_,prob.l);
+ cache = new Cache(prob.l,(long int)(param.cache_size*(1<<20)));
+ QD = new double[prob.l];
+ for(int i=0;i<prob.l;i++)
+ /*** Begin OTB modification ***/
+ QD[i] = (this->*kernel_function)(i,i,m_param);
+ /*** End OTB modification ***/
+ }
+
+ Qfloat *get_Q(int i, int len) const
+ {
+ Qfloat *data;
+ int start, j;
+ if((start = cache->get_data(i,&data,len)) < len)
+ {
+ for(j=start;j<len;j++)
+ /*** Begin OTB modification ***/
+ data[j] = (Qfloat)(y[i]*y[j]*(this->*kernel_function)(i,j,m_param));
+ /*** End OTB modification ***/
+ }
+ return data;
+ }
+
+ double *get_QD() const
+ {
+ return QD;
+ }
+
+ void swap_index(int i, int j) const
+ {
+ cache->swap_index(i,j);
+ Kernel::swap_index(i,j);
+ swap(y[i],y[j]);
+ swap(QD[i],QD[j]);
+ }
+
+ ~SVC_Q()
+ {
+ delete[] y;
+ delete cache;
+ delete[] QD;
+ }
+private:
+ schar *y;
+ Cache *cache;
+ double *QD;
+};
+
+class ONE_CLASS_Q: public Kernel
+{
+public:
+ ONE_CLASS_Q(const svm_problem& prob, const svm_parameter& param)
+ :Kernel(prob.l, prob.x, param)
+ {
+ cache = new Cache(prob.l,(long int)(param.cache_size*(1<<20)));
+ QD = new double[prob.l];
+ for(int i=0;i<prob.l;i++)
+ /*** Begin OTB modification ***/
+ QD[i] = (this->*kernel_function)(i,i,m_param);
+ /*** End OTB modification ***/
+ }
+
+ Qfloat *get_Q(int i, int len) const
+ {
+ Qfloat *data;
+ int start, j;
+ if((start = cache->get_data(i,&data,len)) < len)
+ {
+ for(j=start;j<len;j++)
+ /*** Begin OTB modification ***/
+ data[j] = (Qfloat)(this->*kernel_function)(i,j,m_param);
+ /*** End OTB modification ***/
+ }
+ return data;
+ }
+
+ double *get_QD() const
+ {
+ return QD;
+ }
+
+ void swap_index(int i, int j) const
+ {
+ cache->swap_index(i,j);
+ Kernel::swap_index(i,j);
+ swap(QD[i],QD[j]);
+ }
+
+ ~ONE_CLASS_Q()
+ {
+ delete cache;
+ delete[] QD;
+ }
+private:
+ Cache *cache;
+ double *QD;
+};
+
+class SVR_Q: public Kernel
+{
+public:
+ SVR_Q(const svm_problem& prob, const svm_parameter& param)
+ :Kernel(prob.l, prob.x, param)
+ {
+ l = prob.l;
+ cache = new Cache(l,(long int)(param.cache_size*(1<<20)));
+ QD = new double[2*l];
+ sign = new schar[2*l];
+ index = new int[2*l];
+ for(int k=0;k<l;k++)
+ {
+ sign[k] = 1;
+ sign[k+l] = -1;
+ index[k] = k;
+ index[k+l] = k;
+ /*** Begin OTB modification ***/
+ QD[k] = (this->*kernel_function)(k,k,m_param);
+ /*** End OTB modification ***/
+ QD[k+l] = QD[k];
+ }
+ buffer[0] = new Qfloat[2*l];
+ buffer[1] = new Qfloat[2*l];
+ next_buffer = 0;
+ }
+
+ void swap_index(int i, int j) const
+ {
+ swap(sign[i],sign[j]);
+ swap(index[i],index[j]);
+ swap(QD[i],QD[j]);
+ }
+
+ Qfloat *get_Q(int i, int len) const
+ {
+ Qfloat *data;
+ int j, real_i = index[i];
+ if(cache->get_data(real_i,&data,l) < l)
+ {
+ for(j=0;j<l;j++)
+ /*** Begin OTB modification ***/
+ data[j] = (Qfloat)(this->*kernel_function)(real_i,j,m_param);
+ /*** End OTB modification ***/
+ }
+
+ // reorder and copy
+ Qfloat *buf = buffer[next_buffer];
+ next_buffer = 1 - next_buffer;
+ schar si = sign[i];
+ for(j=0;j<len;j++)
+ buf[j] = (Qfloat) si * (Qfloat) sign[j] * data[index[j]];
+ return buf;
+ }
+
+ double *get_QD() const
+ {
+ return QD;
+ }
+
+ ~SVR_Q()
+ {
+ delete cache;
+ delete[] sign;
+ delete[] index;
+ delete[] buffer[0];
+ delete[] buffer[1];
+ delete[] QD;
+ }
+private:
+ int l;
+ Cache *cache;
+ schar *sign;
+ int *index;
+ mutable int next_buffer;
+ Qfloat *buffer[2];
+ double *QD;
+};
+
+//
+// construct and solve various formulations
+//
+static void solve_c_svc(
+ const svm_problem *prob, const svm_parameter* param,
+ double *alpha, Solver::SolutionInfo* si, double Cp, double Cn)
+{
+ int l = prob->l;
+ double *minus_ones = new double[l];
+ schar *y = new schar[l];
+
+ int i;
+
+ for(i=0;i<l;i++)
+ {
+ alpha[i] = 0;
+ minus_ones[i] = -1;
+ if(prob->y[i] > 0) y[i] = +1; else y[i] = -1;
+ }
+
+ Solver s;
+ s.Solve(l, SVC_Q(*prob,*param,y), minus_ones, y,
+ alpha, Cp, Cn, param->eps, si, param->shrinking);
+
+ double sum_alpha=0;
+ for(i=0;i<l;i++)
+ sum_alpha += alpha[i];
+
+ if (Cp==Cn)
+ info("nu = %f\n", sum_alpha/(Cp*prob->l));
+
+ for(i=0;i<l;i++)
+ alpha[i] *= y[i];
+
+ delete[] minus_ones;
+ delete[] y;
+}
+
+static void solve_nu_svc(
+ const svm_problem *prob, const svm_parameter *param,
+ double *alpha, Solver::SolutionInfo* si)
+{
+ int i;
+ int l = prob->l;
+ double nu = param->nu;
+
+ schar *y = new schar[l];
+
+ for(i=0;i<l;i++)
+ if(prob->y[i]>0)
+ y[i] = +1;
+ else
+ y[i] = -1;
+
+ double sum_pos = nu*l/2;
+ double sum_neg = nu*l/2;
+
+ for(i=0;i<l;i++)
+ if(y[i] == +1)
+ {
+ alpha[i] = min(1.0,sum_pos);
+ sum_pos -= alpha[i];
+ }
+ else
+ {
+ alpha[i] = min(1.0,sum_neg);
+ sum_neg -= alpha[i];
+ }
+
+ double *zeros = new double[l];
+
+ for(i=0;i<l;i++)
+ zeros[i] = 0;
+
+ Solver_NU s;
+ s.Solve(l, SVC_Q(*prob,*param,y), zeros, y,
+ alpha, 1.0, 1.0, param->eps, si, param->shrinking);
+ double r = si->r;
+
+ info("C = %f\n",1/r);
+
+ for(i=0;i<l;i++)
+ alpha[i] *= y[i]/r;
+
+ si->rho /= r;
+ si->obj /= (r*r);
+ si->upper_bound_p = 1/r;
+ si->upper_bound_n = 1/r;
+
+ delete[] y;
+ delete[] zeros;
+}
+
+static void solve_one_class(
+ const svm_problem *prob, const svm_parameter *param,
+ double *alpha, Solver::SolutionInfo* si)
+{
+ int l = prob->l;
+ double *zeros = new double[l];
+ schar *ones = new schar[l];
+ int i;
+
+ int n = (int)(param->nu*prob->l); // # of alpha's at upper bound
+
+ for(i=0;i<n;i++)
+ alpha[i] = 1;
+ if(n<prob->l)
+ alpha[n] = param->nu * prob->l - n;
+ for(i=n+1;i<l;i++)
+ alpha[i] = 0;
+
+ for(i=0;i<l;i++)
+ {
+ zeros[i] = 0;
+ ones[i] = 1;
+ }
+
+ Solver s;
+ s.Solve(l, ONE_CLASS_Q(*prob,*param), zeros, ones,
+ alpha, 1.0, 1.0, param->eps, si, param->shrinking);
+
+ delete[] zeros;
+ delete[] ones;
+}
+
+static void solve_epsilon_svr(
+ const svm_problem *prob, const svm_parameter *param,
+ double *alpha, Solver::SolutionInfo* si)
+{
+ int l = prob->l;
+ double *alpha2 = new double[2*l];
+ double *linear_term = new double[2*l];
+ schar *y = new schar[2*l];
+ int i;
+
+ for(i=0;i<l;i++)
+ {
+ alpha2[i] = 0;
+ linear_term[i] = param->p - prob->y[i];
+ y[i] = 1;
+
+ alpha2[i+l] = 0;
+ linear_term[i+l] = param->p + prob->y[i];
+ y[i+l] = -1;
+ }
+
+ Solver s;
+ s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
+ alpha2, param->C, param->C, param->eps, si, param->shrinking);
+
+ double sum_alpha = 0;
+ for(i=0;i<l;i++)
+ {
+ alpha[i] = alpha2[i] - alpha2[i+l];
+ sum_alpha += fabs(alpha[i]);
+ }
+ info("nu = %f\n",sum_alpha/(param->C*l));
+
+ delete[] alpha2;
+ delete[] linear_term;
+ delete[] y;
+}
+
+static void solve_nu_svr(
+ const svm_problem *prob, const svm_parameter *param,
+ double *alpha, Solver::SolutionInfo* si)
+{
+ int l = prob->l;
+ double C = param->C;
+ double *alpha2 = new double[2*l];
+ double *linear_term = new double[2*l];
+ schar *y = new schar[2*l];
+ int i;
+
+ double sum = C * param->nu * l / 2;
+ for(i=0;i<l;i++)
+ {
+ alpha2[i] = alpha2[i+l] = min(sum,C);
+ sum -= alpha2[i];
+
+ linear_term[i] = - prob->y[i];
+ y[i] = 1;
+
+ linear_term[i+l] = prob->y[i];
+ y[i+l] = -1;
+ }
+
+ Solver_NU s;
+ s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
+ alpha2, C, C, param->eps, si, param->shrinking);
+
+ info("epsilon = %f\n",-si->r);
+
+ for(i=0;i<l;i++)
+ alpha[i] = alpha2[i] - alpha2[i+l];
+
+ delete[] alpha2;
+ delete[] linear_term;
+ delete[] y;
+}
+
+//
+// decision_function
+//
+struct decision_function
+{
+ double *alpha;
+ double rho;
+};
+
+static decision_function svm_train_one(
+ const svm_problem *prob, const svm_parameter *param,
+ double Cp, double Cn)
+{
+ double *alpha = Malloc(double,prob->l);
+ Solver::SolutionInfo si;
+ switch(param->svm_type)
+ {
+ case C_SVC:
+ solve_c_svc(prob,param,alpha,&si,Cp,Cn);
+ break;
+ case NU_SVC:
+ solve_nu_svc(prob,param,alpha,&si);
+ break;
+ case ONE_CLASS:
+ solve_one_class(prob,param,alpha,&si);
+ break;
+ case EPSILON_SVR:
+ solve_epsilon_svr(prob,param,alpha,&si);
+ break;
+ case NU_SVR:
+ solve_nu_svr(prob,param,alpha,&si);
+ break;
+ }
+
+ info("obj = %f, rho = %f\n",si.obj,si.rho);
+
+ // output SVs
+
+ int nSV = 0;
+ int nBSV = 0;
+ for(int i=0;i<prob->l;i++)
+ {
+ if(fabs(alpha[i]) > 0)
+ {
+ ++nSV;
+ if(prob->y[i] > 0)
+ {
+ if(fabs(alpha[i]) >= si.upper_bound_p)
+ ++nBSV;
+ }
+ else
+ {
+ if(fabs(alpha[i]) >= si.upper_bound_n)
+ ++nBSV;
+ }
+ }
+ }
+
+ info("nSV = %d, nBSV = %d\n",nSV,nBSV);
+
+ decision_function f;
+ f.alpha = alpha;
+ f.rho = si.rho;
+ return f;
+}
+
+// Platt's binary SVM Probablistic Output: an improvement from Lin et al.
+static void sigmoid_train(
+ int l, const double *dec_values, const double *labels,
+ double& A, double& B)
+{
+ double prior1=0, prior0 = 0;
+ int i;
+
+ for (i=0;i<l;i++)
+ if (labels[i] > 0) prior1+=1;
+ else prior0+=1;
+
+ int max_iter=100; // Maximal number of iterations
+ double min_step=1e-10; // Minimal step taken in line search
+ double sigma=1e-12; // For numerically strict PD of Hessian
+ double eps=1e-5;
+ double hiTarget=(prior1+1.0)/(prior1+2.0);
+ double loTarget=1/(prior0+2.0);
+ double *t=Malloc(double,l);
+ double fApB,p,q,h11,h22,h21,g1,g2,det,dA,dB,gd,stepsize;
+ double newA,newB,newf,d1,d2;
+ int iter;
+
+ // Initial Point and Initial Fun Value
+ A=0.0; B=log((prior0+1.0)/(prior1+1.0));
+ double fval = 0.0;
+
+ for (i=0;i<l;i++)
+ {
+ if (labels[i]>0) t[i]=hiTarget;
+ else t[i]=loTarget;
+ fApB = dec_values[i]*A+B;
+ if (fApB>=0)
+ fval += t[i]*fApB + log(1+exp(-fApB));
+ else
+ fval += (t[i] - 1)*fApB +log(1+exp(fApB));
+ }
+ for (iter=0;iter<max_iter;iter++)
+ {
+ // Update Gradient and Hessian (use H' = H + sigma I)
+ h11=sigma; // numerically ensures strict PD
+ h22=sigma;
+ h21=0.0;g1=0.0;g2=0.0;
+ for (i=0;i<l;i++)
+ {
+ fApB = dec_values[i]*A+B;
+ if (fApB >= 0)
+ {
+ p=exp(-fApB)/(1.0+exp(-fApB));
+ q=1.0/(1.0+exp(-fApB));
+ }
+ else
+ {
+ p=1.0/(1.0+exp(fApB));
+ q=exp(fApB)/(1.0+exp(fApB));
+ }
+ d2=p*q;
+ h11+=dec_values[i]*dec_values[i]*d2;
+ h22+=d2;
+ h21+=dec_values[i]*d2;
+ d1=t[i]-p;
+ g1+=dec_values[i]*d1;
+ g2+=d1;
+ }
+
+ // Stopping Criteria
+ if (fabs(g1)<eps && fabs(g2)<eps)
+ break;
+
+ // Finding Newton direction: -inv(H') * g
+ det=h11*h22-h21*h21;
+ dA=-(h22*g1 - h21 * g2) / det;
+ dB=-(-h21*g1+ h11 * g2) / det;
+ gd=g1*dA+g2*dB;
+
+
+ stepsize = 1; // Line Search
+ while (stepsize >= min_step)
+ {
+ newA = A + stepsize * dA;
+ newB = B + stepsize * dB;
+
+ // New function value
+ newf = 0.0;
+ for (i=0;i<l;i++)
+ {
+ fApB = dec_values[i]*newA+newB;
+ if (fApB >= 0)
+ newf += t[i]*fApB + log(1+exp(-fApB));
+ else
+ newf += (t[i] - 1)*fApB +log(1+exp(fApB));
+ }
+ // Check sufficient decrease
+ if (newf<fval+0.0001*stepsize*gd)
+ {
+ A=newA;B=newB;fval=newf;
+ break;
+ }
+ else
+ stepsize = stepsize / 2.0;
+ }
+
+ if (stepsize < min_step)
+ {
+ info("Line search fails in two-class probability estimates\n");
+ break;
+ }
+ }
+
+ if (iter>=max_iter)
+ info("Reaching maximal iterations in two-class probability estimates\n");
+ free(t);
+}
+
+static double sigmoid_predict(double decision_value, double A, double B)
+{
+ double fApB = decision_value*A+B;
+ if (fApB >= 0)
+ return exp(-fApB)/(1.0+exp(-fApB));
+ else
+ return 1.0/(1+exp(fApB)) ;
+}
+
+// Method 2 from the multiclass_prob paper by Wu, Lin, and Weng
+static void multiclass_probability(int k, double **r, double *p)
+{
+ int t,j;
+ int iter = 0, max_iter=max(100,k);
+ double **Q=Malloc(double *,k);
+ double *Qp=Malloc(double,k);
+ double pQp, eps=0.005/k;
+
+ for (t=0;t<k;t++)
+ {
+ p[t]=1.0/k; // Valid if k = 1
+ Q[t]=Malloc(double,k);
+ Q[t][t]=0;
+ for (j=0;j<t;j++)
+ {
+ Q[t][t]+=r[j][t]*r[j][t];
+ Q[t][j]=Q[j][t];
+ }
+ for (j=t+1;j<k;j++)
+ {
+ Q[t][t]+=r[j][t]*r[j][t];
+ Q[t][j]=-r[j][t]*r[t][j];
+ }
+ }
+ for (iter=0;iter<max_iter;iter++)
+ {
+ // stopping condition, recalculate QP,pQP for numerical accuracy
+ pQp=0;
+ for (t=0;t<k;t++)
+ {
+ Qp[t]=0;
+ for (j=0;j<k;j++)
+ Qp[t]+=Q[t][j]*p[j];
+ pQp+=p[t]*Qp[t];
+ }
+ double max_error=0;
+ for (t=0;t<k;t++)
+ {
+ double error=fabs(Qp[t]-pQp);
+ if (error>max_error)
+ max_error=error;
+ }
+ if (max_error<eps) break;
+
+ for (t=0;t<k;t++)
+ {
+ double diff=(-Qp[t]+pQp)/Q[t][t];
+ p[t]+=diff;
+ pQp=(pQp+diff*(diff*Q[t][t]+2*Qp[t]))/(1+diff)/(1+diff);
+ for (j=0;j<k;j++)
+ {
+ Qp[j]=(Qp[j]+diff*Q[t][j])/(1+diff);
+ p[j]/=(1+diff);
+ }
+ }
+ }
+ if (iter>=max_iter)
+ info("Exceeds max_iter in multiclass_prob\n");
+ for(t=0;t<k;t++) free(Q[t]);
+ free(Q);
+ free(Qp);
+}
+
+// Cross-validation decision values for probability estimates
+static void svm_binary_svc_probability(
+ const svm_problem *prob, const svm_parameter *param,
+ double Cp, double Cn, double& probA, double& probB)
+{
+ int i;
+ int nr_fold = 5;
+ int *perm = Malloc(int,prob->l);
+ double *dec_values = Malloc(double,prob->l);
+
+ // random shuffle
+ for(i=0;i<prob->l;i++) perm[i]=i;
+ for(i=0;i<prob->l;i++)
+ {
+ int j = i+rand()%(prob->l-i);
+ swap(perm[i],perm[j]);
+ }
+ for(i=0;i<nr_fold;i++)
+ {
+ int begin = i*prob->l/nr_fold;
+ int end = (i+1)*prob->l/nr_fold;
+ int j,k;
+ struct svm_problem subprob;
+
+ subprob.l = prob->l-(end-begin);
+ subprob.x = Malloc(struct svm_node*,subprob.l);
+ subprob.y = Malloc(double,subprob.l);
+
+ k=0;
+ for(j=0;j<begin;j++)
+ {
+ subprob.x[k] = prob->x[perm[j]];
+ subprob.y[k] = prob->y[perm[j]];
+ ++k;
+ }
+ for(j=end;j<prob->l;j++)
+ {
+ subprob.x[k] = prob->x[perm[j]];
+ subprob.y[k] = prob->y[perm[j]];
+ ++k;
+ }
+ int p_count=0,n_count=0;
+ for(j=0;j<k;j++)
+ if(subprob.y[j]>0)
+ p_count++;
+ else
+ n_count++;
+
+ if(p_count==0 && n_count==0)
+ for(j=begin;j<end;j++)
+ dec_values[perm[j]] = 0;
+ else if(p_count > 0 && n_count == 0)
+ for(j=begin;j<end;j++)
+ dec_values[perm[j]] = 1;
+ else if(p_count == 0 && n_count > 0)
+ for(j=begin;j<end;j++)
+ dec_values[perm[j]] = -1;
+ else
+ {
+ svm_parameter subparam = *param;
+ subparam.probability=0;
+ subparam.C=1.0;
+ subparam.nr_weight=2;
+ subparam.weight_label = Malloc(int,2);
+ subparam.weight = Malloc(double,2);
+ subparam.weight_label[0]=+1;
+ subparam.weight_label[1]=-1;
+ subparam.weight[0]=Cp;
+ subparam.weight[1]=Cn;
+ struct svm_model *submodel = svm_train(&subprob,&subparam);
+ for(j=begin;j<end;j++)
+ {
+ svm_predict_values(submodel,prob->x[perm[j]],&(dec_values[perm[j]]));
+ // ensure +1 -1 order; reason not using CV subroutine
+ dec_values[perm[j]] *= submodel->label[0];
+ }
+ svm_free_and_destroy_model(&submodel);
+ svm_destroy_param(&subparam);
+ }
+ free(subprob.x);
+ free(subprob.y);
+ }
+ sigmoid_train(prob->l,dec_values,prob->y,probA,probB);
+ free(dec_values);
+ free(perm);
+}
+
+// Return parameter of a Laplace distribution
+static double svm_svr_probability(
+ const svm_problem *prob, const svm_parameter *param)
+{
+ int i;
+ int nr_fold = 5;
+ double *ymv = Malloc(double,prob->l);
+ double mae = 0;
+
+ svm_parameter newparam = *param;
+ newparam.probability = 0;
+ svm_cross_validation(prob,&newparam,nr_fold,ymv);
+ for(i=0;i<prob->l;i++)
+ {
+ ymv[i]=prob->y[i]-ymv[i];
+ mae += fabs(ymv[i]);
+ }
+ mae /= prob->l;
+ double std=sqrt(2*mae*mae);
+ int count=0;
+ mae=0;
+ for(i=0;i<prob->l;i++)
+ if (fabs(ymv[i]) > 5*std)
+ count=count+1;
+ else
+ mae+=fabs(ymv[i]);
+ mae /= (prob->l-count);
+ info("Prob. model for test data: target value = predicted value + z,\nz: Laplace distribution e^(-|z|/sigma)/(2sigma),sigma= %g\n",mae);
+ free(ymv);
+ return mae;
+}
+
+
+// label: label name, start: begin of each class, count: #data of classes, perm: indices to the original data
+// perm, length l, must be allocated before calling this subroutine
+static void svm_group_classes(const svm_problem *prob, int *nr_class_ret, int **label_ret, int **start_ret, int **count_ret, int *perm)
+{
+ int l = prob->l;
+ int max_nr_class = 16;
+ int nr_class = 0;
+ int *label = Malloc(int,max_nr_class);
+ int *count = Malloc(int,max_nr_class);
+ int *data_label = Malloc(int,l);
+ int i;
+
+ for(i=0;i<l;i++)
+ {
+ int this_label = (int)prob->y[i];
+ int j;
+ for(j=0;j<nr_class;j++)
+ {
+ if(this_label == label[j])
+ {
+ ++count[j];
+ break;
+ }
+ }
+ data_label[i] = j;
+ if(j == nr_class)
+ {
+ if(nr_class == max_nr_class)
+ {
+ max_nr_class *= 2;
+ label = (int *)realloc(label,max_nr_class*sizeof(int));
+ count = (int *)realloc(count,max_nr_class*sizeof(int));
+ }
+ label[nr_class] = this_label;
+ count[nr_class] = 1;
+ ++nr_class;
+ }
+ }
+
+ int *start = Malloc(int,nr_class);
+ start[0] = 0;
+ for(i=1;i<nr_class;i++)
+ start[i] = start[i-1]+count[i-1];
+ for(i=0;i<l;i++)
+ {
+ perm[start[data_label[i]]] = i;
+ ++start[data_label[i]];
+ }
+ start[0] = 0;
+ for(i=1;i<nr_class;i++)
+ start[i] = start[i-1]+count[i-1];
+
+ *nr_class_ret = nr_class;
+ *label_ret = label;
+ *start_ret = start;
+ *count_ret = count;
+ free(data_label);
+}
+
+//
+// Interface functions
+//
+svm_model *svm_train(const svm_problem *prob, const svm_parameter *param)
+{
+ svm_model *model = Malloc(svm_model,1);
+ model->param = *param;
+ model->free_sv = 0; // XXX
+
+ if(param->svm_type == ONE_CLASS ||
+ param->svm_type == EPSILON_SVR ||
+ param->svm_type == NU_SVR)
+ {
+ // regression or one-class-svm
+ model->nr_class = 2;
+ model->label = NULL;
+ model->nSV = NULL;
+ model->probA = NULL; model->probB = NULL;
+ model->sv_coef = Malloc(double *,1);
+
+ if(param->probability &&
+ (param->svm_type == EPSILON_SVR ||
+ param->svm_type == NU_SVR))
+ {
+ model->probA = Malloc(double,1);
+ model->probA[0] = svm_svr_probability(prob,param);
+ }
+
+ decision_function f = svm_train_one(prob,param,0,0);
+ model->rho = Malloc(double,1);
+ model->rho[0] = f.rho;
+
+ int nSV = 0;
+ int i;
+ for(i=0;i<prob->l;i++)
+ if(fabs(f.alpha[i]) > 0) ++nSV;
+ model->l = nSV;
+ model->SV = Malloc(svm_node *,nSV);
+ model->sv_coef[0] = Malloc(double,nSV);
+ int j = 0;
+ for(i=0;i<prob->l;i++)
+ if(fabs(f.alpha[i]) > 0)
+ {
+ model->SV[j] = prob->x[i];
+ model->sv_coef[0][j] = f.alpha[i];
+ ++j;
+ }
+
+ free(f.alpha);
+ }
+ else
+ {
+ // classification
+ int l = prob->l;
+ int nr_class;
+ int *label = NULL;
+ int *start = NULL;
+ int *count = NULL;
+ int *perm = Malloc(int,l);
+
+ // group training data of the same class
+ svm_group_classes(prob,&nr_class,&label,&start,&count,perm);
+ svm_node **x = Malloc(svm_node *,l);
+ int i;
+ for(i=0;i<l;i++)
+ x[i] = prob->x[perm[i]];
+
+ // calculate weighted C
+
+ double *weighted_C = Malloc(double, nr_class);
+ for(i=0;i<nr_class;i++)
+ weighted_C[i] = param->C;
+ for(i=0;i<param->nr_weight;i++)
+ {
+ int j;
+ for(j=0;j<nr_class;j++)
+ if(param->weight_label[i] == label[j])
+ break;
+ if(j == nr_class)
+ fprintf(stderr,"warning: class label %d specified in weight is not found\n", param->weight_label[i]);
+ else
+ weighted_C[j] *= param->weight[i];
+ }
+
+ // train k*(k-1)/2 models
+
+ bool *nonzero = Malloc(bool,l);
+ for(i=0;i<l;i++)
+ nonzero[i] = false;
+ decision_function *f = Malloc(decision_function,nr_class*(nr_class-1)/2);
+
+ double *probA=NULL,*probB=NULL;
+ if (param->probability)
+ {
+ probA=Malloc(double,nr_class*(nr_class-1)/2);
+ probB=Malloc(double,nr_class*(nr_class-1)/2);
+ }
+
+ int p = 0;
+ for(i=0;i<nr_class;i++)
+ for(int j=i+1;j<nr_class;j++)
+ {
+ svm_problem sub_prob;
+ int si = start[i], sj = start[j];
+ int ci = count[i], cj = count[j];
+ sub_prob.l = ci+cj;
+ sub_prob.x = Malloc(svm_node *,sub_prob.l);
+ sub_prob.y = Malloc(double,sub_prob.l);
+ int k;
+ for(k=0;k<ci;k++)
+ {
+ sub_prob.x[k] = x[si+k];
+ sub_prob.y[k] = +1;
+ }
+ for(k=0;k<cj;k++)
+ {
+ sub_prob.x[ci+k] = x[sj+k];
+ sub_prob.y[ci+k] = -1;
+ }
+
+ if(param->probability)
+ svm_binary_svc_probability(&sub_prob,param,weighted_C[i],weighted_C[j],probA[p],probB[p]);
+
+ f[p] = svm_train_one(&sub_prob,param,weighted_C[i],weighted_C[j]);
+ for(k=0;k<ci;k++)
+ if(!nonzero[si+k] && fabs(f[p].alpha[k]) > 0)
+ nonzero[si+k] = true;
+ for(k=0;k<cj;k++)
+ if(!nonzero[sj+k] && fabs(f[p].alpha[ci+k]) > 0)
+ nonzero[sj+k] = true;
+ free(sub_prob.x);
+ free(sub_prob.y);
+ ++p;
+ }
+
+ // build output
+
+ model->nr_class = nr_class;
+
+ model->label = Malloc(int,nr_class);
+ for(i=0;i<nr_class;i++)
+ model->label[i] = label[i];
+
+ model->rho = Malloc(double,nr_class*(nr_class-1)/2);
+ for(i=0;i<nr_class*(nr_class-1)/2;i++)
+ model->rho[i] = f[i].rho;
+
+ if(param->probability)
+ {
+ model->probA = Malloc(double,nr_class*(nr_class-1)/2);
+ model->probB = Malloc(double,nr_class*(nr_class-1)/2);
+ for(i=0;i<nr_class*(nr_class-1)/2;i++)
+ {
+ model->probA[i] = probA[i];
+ model->probB[i] = probB[i];
+ }
+ }
+ else
+ {
+ model->probA=NULL;
+ model->probB=NULL;
+ }
+
+ int total_sv = 0;
+ int *nz_count = Malloc(int,nr_class);
+ model->nSV = Malloc(int,nr_class);
+ for(i=0;i<nr_class;i++)
+ {
+ int nSV = 0;
+ for(int j=0;j<count[i];j++)
+ if(nonzero[start[i]+j])
+ {
+ ++nSV;
+ ++total_sv;
+ }
+ model->nSV[i] = nSV;
+ nz_count[i] = nSV;
+ }
+
+ info("Total nSV = %d\n",total_sv);
+
+ model->l = total_sv;
+ model->SV = Malloc(svm_node *,total_sv);
+ p = 0;
+ for(i=0;i<l;i++)
+ if(nonzero[i]) model->SV[p++] = x[i];
+
+ int *nz_start = Malloc(int,nr_class);
+ nz_start[0] = 0;
+ for(i=1;i<nr_class;i++)
+ nz_start[i] = nz_start[i-1]+nz_count[i-1];
+
+ model->sv_coef = Malloc(double *,nr_class-1);
+ for(i=0;i<nr_class-1;i++)
+ model->sv_coef[i] = Malloc(double,total_sv);
+
+ p = 0;
+ for(i=0;i<nr_class;i++)
+ for(int j=i+1;j<nr_class;j++)
+ {
+ // classifier (i,j): coefficients with
+ // i are in sv_coef[j-1][nz_start[i]...],
+ // j are in sv_coef[i][nz_start[j]...]
+
+ int si = start[i];
+ int sj = start[j];
+ int ci = count[i];
+ int cj = count[j];
+
+ int q = nz_start[i];
+ int k;
+ for(k=0;k<ci;k++)
+ if(nonzero[si+k])
+ model->sv_coef[j-1][q++] = f[p].alpha[k];
+ q = nz_start[j];
+ for(k=0;k<cj;k++)
+ if(nonzero[sj+k])
+ model->sv_coef[i][q++] = f[p].alpha[ci+k];
+ ++p;
+ }
+
+ free(label);
+ free(probA);
+ free(probB);
+ free(count);
+ free(perm);
+ free(start);
+ free(x);
+ free(weighted_C);
+ free(nonzero);
+ for(i=0;i<nr_class*(nr_class-1)/2;i++)
+ free(f[i].alpha);
+ free(f);
+ free(nz_count);
+ free(nz_start);
+ }
+ return model;
+}
+
+// Stratified cross validation
+void svm_cross_validation(const svm_problem *prob, const svm_parameter *param, int nr_fold, double *target)
+{
+ int i;
+ int *fold_start = Malloc(int,nr_fold+1);
+ int l = prob->l;
+ int *perm = Malloc(int,l);
+ int nr_class;
+
+ // stratified cv may not give leave-one-out rate
+ // Each class to l folds -> some folds may have zero elements
+ if((param->svm_type == C_SVC ||
+ param->svm_type == NU_SVC) && nr_fold < l)
+ {
+ int *start = NULL;
+ int *label = NULL;
+ int *count = NULL;
+ svm_group_classes(prob,&nr_class,&label,&start,&count,perm);
+
+ // random shuffle and then data grouped by fold using the array perm
+ int *fold_count = Malloc(int,nr_fold);
+ int c;
+ int *index = Malloc(int,l);
+ for(i=0;i<l;i++)
+ index[i]=perm[i];
+ for (c=0; c<nr_class; c++)
+ for(i=0;i<count[c];i++)
+ {
+ int j = i+rand()%(count[c]-i);
+ swap(index[start[c]+j],index[start[c]+i]);
+ }
+ for(i=0;i<nr_fold;i++)
+ {
+ fold_count[i] = 0;
+ for (c=0; c<nr_class;c++)
+ fold_count[i]+=(i+1)*count[c]/nr_fold-i*count[c]/nr_fold;
+ }
+ fold_start[0]=0;
+ for (i=1;i<=nr_fold;i++)
+ fold_start[i] = fold_start[i-1]+fold_count[i-1];
+ for (c=0; c<nr_class;c++)
+ for(i=0;i<nr_fold;i++)
+ {
+ int begin = start[c]+i*count[c]/nr_fold;
+ int end = start[c]+(i+1)*count[c]/nr_fold;
+ for(int j=begin;j<end;j++)
+ {
+ perm[fold_start[i]] = index[j];
+ fold_start[i]++;
+ }
+ }
+ fold_start[0]=0;
+ for (i=1;i<=nr_fold;i++)
+ fold_start[i] = fold_start[i-1]+fold_count[i-1];
+ free(start);
+ free(label);
+ free(count);
+ free(index);
+ free(fold_count);
+ }
+ else
+ {
+ for(i=0;i<l;i++) perm[i]=i;
+ for(i=0;i<l;i++)
+ {
+ int j = i+rand()%(l-i);
+ swap(perm[i],perm[j]);
+ }
+ for(i=0;i<=nr_fold;i++)
+ fold_start[i]=i*l/nr_fold;
+ }
+
+ for(i=0;i<nr_fold;i++)
+ {
+ int begin = fold_start[i];
+ int end = fold_start[i+1];
+ int j,k;
+ struct svm_problem subprob;
+
+ subprob.l = l-(end-begin);
+ subprob.x = Malloc(struct svm_node*,subprob.l);
+ subprob.y = Malloc(double,subprob.l);
+
+ k=0;
+ for(j=0;j<begin;j++)
+ {
+ subprob.x[k] = prob->x[perm[j]];
+ subprob.y[k] = prob->y[perm[j]];
+ ++k;
+ }
+ for(j=end;j<l;j++)
+ {
+ subprob.x[k] = prob->x[perm[j]];
+ subprob.y[k] = prob->y[perm[j]];
+ ++k;
+ }
+ struct svm_model *submodel = svm_train(&subprob,param);
+ if(param->probability &&
+ (param->svm_type == C_SVC || param->svm_type == NU_SVC))
+ {
+ double *prob_estimates=Malloc(double,svm_get_nr_class(submodel));
+ for(j=begin;j<end;j++)
+ target[perm[j]] = svm_predict_probability(submodel,prob->x[perm[j]],prob_estimates);
+ free(prob_estimates);
+ }
+ else
+ for(j=begin;j<end;j++)
+ target[perm[j]] = svm_predict(submodel,prob->x[perm[j]]);
+ svm_free_and_destroy_model(&submodel);
+ free(subprob.x);
+ free(subprob.y);
+ }
+ free(fold_start);
+ free(perm);
+}
+
+
+int svm_get_svm_type(const svm_model *model)
+{
+ return model->param.svm_type;
+}
+
+int svm_get_nr_class(const svm_model *model)
+{
+ return model->nr_class;
+}
+
+void svm_get_labels(const svm_model *model, int* label)
+{
+ if (model->label != NULL)
+ for(int i=0;i<model->nr_class;i++)
+ label[i] = model->label[i];
+}
+
+double svm_get_svr_probability(const svm_model *model)
+{
+ if ((model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) &&
+ model->probA!=NULL)
+ return model->probA[0];
+ else
+ {
+ fprintf(stderr,"Model doesn't contain information for SVR probability inference\n");
+ return 0;
+ }
+}
+
+double svm_predict_values(const svm_model *model, const svm_node *x, double* dec_values)
+{
+ if(model->param.svm_type == ONE_CLASS ||
+ model->param.svm_type == EPSILON_SVR ||
+ model->param.svm_type == NU_SVR)
+ {
+ double *sv_coef = model->sv_coef[0];
+ double sum = 0;
+ for(int i=0;i<model->l;i++)
+ sum += sv_coef[i] * Kernel::k_function(x,model->SV[i],model->param);
+ sum -= model->rho[0];
+ *dec_values = sum;
+
+ if(model->param.svm_type == ONE_CLASS)
+ return (sum>0)?1:-1;
+ else
+ return sum;
+ }
+ else
+ {
+ int i;
+ int nr_class = model->nr_class;
+ int l = model->l;
+
+ double *kvalue = Malloc(double,l);
+ for(i=0;i<l;i++)
+ kvalue[i] = Kernel::k_function(x,model->SV[i],model->param);
+
+ int *start = Malloc(int,nr_class);
+ start[0] = 0;
+ for(i=1;i<nr_class;i++)
+ start[i] = start[i-1]+model->nSV[i-1];
+
+ int *vote = Malloc(int,nr_class);
+ for(i=0;i<nr_class;i++)
+ vote[i] = 0;
+
+ int p=0;
+ for(i=0;i<nr_class;i++)
+ for(int j=i+1;j<nr_class;j++)
+ {
+ double sum = 0;
+ int si = start[i];
+ int sj = start[j];
+ int ci = model->nSV[i];
+ int cj = model->nSV[j];
+
+ int k;
+ double *coef1 = model->sv_coef[j-1];
+ double *coef2 = model->sv_coef[i];
+ for(k=0;k<ci;k++)
+ sum += coef1[si+k] * kvalue[si+k];
+ for(k=0;k<cj;k++)
+ sum += coef2[sj+k] * kvalue[sj+k];
+ sum -= model->rho[p];
+ dec_values[p] = sum;
+
+ if(dec_values[p] > 0)
+ ++vote[i];
+ else
+ ++vote[j];
+ p++;
+ }
+
+ int vote_max_idx = 0;
+ for(i=1;i<nr_class;i++)
+ if(vote[i] > vote[vote_max_idx])
+ vote_max_idx = i;
+
+ free(kvalue);
+ free(start);
+ free(vote);
+ return model->label[vote_max_idx];
+ }
+}
+
+double svm_predict(const svm_model *model, const svm_node *x)
+{
+ int nr_class = model->nr_class;
+ double *dec_values;
+ if(model->param.svm_type == ONE_CLASS ||
+ model->param.svm_type == EPSILON_SVR ||
+ model->param.svm_type == NU_SVR)
+ dec_values = Malloc(double, 1);
+ else
+ dec_values = Malloc(double, nr_class*(nr_class-1)/2);
+ double pred_result = svm_predict_values(model, x, dec_values);
+ free(dec_values);
+ return pred_result;
+}
+
+double svm_predict_probability(
+ const svm_model *model, const svm_node *x, double *prob_estimates)
+{
+ if ((model->param.svm_type == C_SVC || model->param.svm_type == NU_SVC) &&
+ model->probA!=NULL && model->probB!=NULL)
+ {
+ int i;
+ int nr_class = model->nr_class;
+ double *dec_values = Malloc(double, nr_class*(nr_class-1)/2);
+ svm_predict_values(model, x, dec_values);
+
+ double min_prob=1e-7;
+ double **pairwise_prob=Malloc(double *,nr_class);
+ for(i=0;i<nr_class;i++)
+ pairwise_prob[i]=Malloc(double,nr_class);
+ int k=0;
+ for(i=0;i<nr_class;i++)
+ for(int j=i+1;j<nr_class;j++)
+ {
+ pairwise_prob[i][j]=min(max(sigmoid_predict(dec_values[k],model->probA[k],model->probB[k]),min_prob),1-min_prob);
+ pairwise_prob[j][i]=1-pairwise_prob[i][j];
+ k++;
+ }
+ multiclass_probability(nr_class,pairwise_prob,prob_estimates);
+
+ int prob_max_idx = 0;
+ for(i=1;i<nr_class;i++)
+ if(prob_estimates[i] > prob_estimates[prob_max_idx])
+ prob_max_idx = i;
+ for(i=0;i<nr_class;i++)
+ free(pairwise_prob[i]);
+ free(dec_values);
+ free(pairwise_prob);
+ return model->label[prob_max_idx];
+ }
+ else
+ return svm_predict(model, x);
+}
+
+static const char *svm_type_table[] =
+{
+ "c_svc","nu_svc","one_class","epsilon_svr","nu_svr",NULL
+};
+
+static const char *kernel_type_table[]=
+{
+ /*** Begin OTB modification ***/
+ "linear","polynomial","rbf","sigmoid","precomputed","generic","composed",NULL
+ /*** End OTB modification ***/
+};
+
+int svm_save_model(const char *model_file_name, const svm_model *model)
+{
+ FILE *fp = fopen(model_file_name,"w");
+ if(fp==NULL) return -1;
+
+ const svm_parameter& param = model->param;
+
+ fprintf(fp,"svm_type %s\n", svm_type_table[param.svm_type]);
+ fprintf(fp,"kernel_type %s\n", kernel_type_table[param.kernel_type]);
+
+ /*** Begin OTB modification ***/
+ if( param.kernel_type == GENERIC )
+ {
+ if( param.kernel_generic == NULL )
+ {
+ fprintf(stderr,"generic kernel functor is not initialized\n");
+ return -1;
+ }
+
+ //Load generic parameters
+ int cr = param.kernel_generic->save_parameters(&fp,"generic_kernel_parameters");
+ if( cr != 0 )
+ {
+ fprintf(stderr,"error while saving generic kernel parameters to the file %s.\n",model_file_name);
+ }
+ }
+ if (param.kernel_type == COMPOSED)
+ {
+ if (param.kernel_composed == NULL)
+ {
+ fprintf(stderr, "composed kernel functor is not initialized\n");
+ return -1;
+ }
+ //Load generic parameters
+ int cr = param.kernel_composed->save_parameters(&fp, "composed_kernel_parameters");
+ if (cr != 0)
+ {
+ fprintf(stderr, "error while saving composed kernel parameters to the file %s.\n", model_file_name);
+ }
+ }
+ /*** End OTB modification ***/
+
+ if(param.kernel_type == POLY)
+ fprintf(fp,"degree %d\n", param.degree);
+
+ if(param.kernel_type == POLY || param.kernel_type == RBF || param.kernel_type == SIGMOID)
+ fprintf(fp,"gamma %g\n", param.gamma);
+
+ if(param.kernel_type == POLY || param.kernel_type == SIGMOID)
+ fprintf(fp,"coef0 %g\n", param.coef0);
+
+ int nr_class = model->nr_class;
+ int l = model->l;
+ fprintf(fp, "nr_class %d\n", nr_class);
+ fprintf(fp, "total_sv %d\n",l);
+
+ {
+ fprintf(fp, "rho");
+ for(int i=0;i<nr_class*(nr_class-1)/2;i++)
+ fprintf(fp," %g",model->rho[i]);
+ fprintf(fp, "\n");
+ }
+
+ if(model->label)
+ {
+ fprintf(fp, "label");
+ for(int i=0;i<nr_class;i++)
+ fprintf(fp," %d",model->label[i]);
+ fprintf(fp, "\n");
+ }
+
+ if(model->probA) // regression has probA only
+ {
+ fprintf(fp, "probA");
+ for(int i=0;i<nr_class*(nr_class-1)/2;i++)
+ fprintf(fp," %g",model->probA[i]);
+ fprintf(fp, "\n");
+ }
+ if(model->probB)
+ {
+ fprintf(fp, "probB");
+ for(int i=0;i<nr_class*(nr_class-1)/2;i++)
+ fprintf(fp," %g",model->probB[i]);
+ fprintf(fp, "\n");
+ }
+
+ if(model->nSV)
+ {
+ fprintf(fp, "nr_sv");
+ for(int i=0;i<nr_class;i++)
+ fprintf(fp," %d",model->nSV[i]);
+ fprintf(fp, "\n");
+ }
+
+ fprintf(fp, "SV\n");
+ const double * const *sv_coef = model->sv_coef;
+ const svm_node * const *SV = model->SV;
+
+ for(int i=0;i<l;i++)
+ {
+ for(int j=0;j<nr_class-1;j++)
+ fprintf(fp, "%.16g ",sv_coef[j][i]);
+
+ const svm_node *p = SV[i];
+
+ if(param.kernel_type == PRECOMPUTED)
+ fprintf(fp,"0:%d ",(int)(p->value));
+ else
+ while(p->index != -1)
+ {
+ fprintf(fp,"%d:%.8g ",p->index,p->value);
+ p++;
+ }
+ fprintf(fp, "\n");
+ }
+ if (ferror(fp) != 0 || fclose(fp) != 0) return -1;
+ else return 0;
+}
+
+static char *line = NULL;
+static int max_line_len;
+
+static char* readline(FILE *input)
+{
+ int len;
+
+ if(fgets(line,max_line_len,input) == NULL)
+ return NULL;
+
+ while(strrchr(line,'\n') == NULL)
+ {
+ max_line_len *= 2;
+ line = (char *) realloc(line,max_line_len);
+ len = (int) strlen(line);
+ if(fgets(line+len,max_line_len-len,input) == NULL)
+ break;
+ }
+ return line;
+}
+
+svm_model *svm_load_model(const char *model_file_name, GenericKernelFunctorBase* generic_kernel_functor)
+{
+ FILE *fp = fopen(model_file_name,"rb");
+ if(fp==NULL) return NULL;
+
+ // read parameters
+
+ svm_model *model = Malloc(svm_model,1);
+ svm_parameter& param = model->param;
+ model->rho = NULL;
+ model->probA = NULL;
+ model->probB = NULL;
+ model->label = NULL;
+ model->nSV = NULL;
+
+ char cmd[81];
+ while(1)
+ {
+ fscanf(fp,"%80s",cmd);
+
+ if(strcmp(cmd,"svm_type")==0)
+ {
+ fscanf(fp,"%80s",cmd);
+ int i;
+ for(i=0;svm_type_table[i];i++)
+ {
+ if(strcmp(svm_type_table[i],cmd)==0)
+ {
+ param.svm_type=i;
+ break;
+ }
+ }
+ if(svm_type_table[i] == NULL)
+ {
+ fprintf(stderr,"unknown svm type.\n");
+ free(model->rho);
+ free(model->label);
+ free(model->nSV);
+ free(model);
+ return NULL;
+ }
+ }
+ else if(strcmp(cmd,"kernel_type")==0)
+ {
+ fscanf(fp,"%80s",cmd);
+ int i;
+ for(i=0;kernel_type_table[i];i++)
+ {
+ if(strcmp(kernel_type_table[i],cmd)==0)
+ {
+ param.kernel_type=i;
+ break;
+ }
+ }
+ if(kernel_type_table[i] == NULL)
+ {
+ fprintf(stderr,"unknown kernel function.\n");
+ free(model->rho);
+ free(model->label);
+ free(model->nSV);
+ free(model);
+ return NULL;
+ }
+ }
+ else if(strcmp(cmd,"degree")==0)
+ fscanf(fp,"%d",¶m.degree);
+ else if(strcmp(cmd,"gamma")==0)
+ fscanf(fp,"%lf",¶m.gamma);
+ else if(strcmp(cmd,"coef0")==0)
+ fscanf(fp,"%lf",¶m.coef0);
+ else if(strcmp(cmd,"nr_class")==0)
+ fscanf(fp,"%d",&model->nr_class);
+ else if(strcmp(cmd,"total_sv")==0)
+ fscanf(fp,"%d",&model->l);
+ else if(strcmp(cmd,"rho")==0)
+ {
+ int n = model->nr_class * (model->nr_class-1)/2;
+ model->rho = Malloc(double,n);
+ for(int i=0;i<n;i++)
+ fscanf(fp,"%lf",&model->rho[i]);
+ }
+ else if(strcmp(cmd,"label")==0)
+ {
+ int n = model->nr_class;
+ model->label = Malloc(int,n);
+ for(int i=0;i<n;i++)
+ fscanf(fp,"%d",&model->label[i]);
+ }
+ else if(strcmp(cmd,"probA")==0)
+ {
+ int n = model->nr_class * (model->nr_class-1)/2;
+ model->probA = Malloc(double,n);
+ for(int i=0;i<n;i++)
+ fscanf(fp,"%lf",&model->probA[i]);
+ }
+ else if(strcmp(cmd,"probB")==0)
+ {
+ int n = model->nr_class * (model->nr_class-1)/2;
+ model->probB = Malloc(double,n);
+ for(int i=0;i<n;i++)
+ fscanf(fp,"%lf",&model->probB[i]);
+ }
+ else if(strcmp(cmd,"nr_sv")==0)
+ {
+ int n = model->nr_class;
+ model->nSV = Malloc(int,n);
+ for(int i=0;i<n;i++)
+ fscanf(fp,"%d",&model->nSV[i]);
+ }
+ else if(strcmp(cmd,"SV")==0)
+ {
+ while(1)
+ {
+ int c = getc(fp);
+ if(c==EOF || c=='\n') break;
+ }
+ break;
+ }
+ /*** Begin OTB modification ***/
+ else if (strcmp(cmd, "generic_kernel_parameters") == 0)
+ {
+ if (param.kernel_type == GENERIC)
+ {
+ if (generic_kernel_functor == NULL)
+ {
+ fprintf(stderr, "generic kernel functor is not initialized\n", cmd);
+ return NULL;
+ }
+ param.kernel_generic = generic_kernel_functor;
+ //Load generic parameters
+ int cr = param.kernel_generic->load_parameters(&fp);
+ if (cr != 0)
+ {
+ fprintf(stderr, "error while loading generic kernel parameters from the file %s.\n", model_file_name);
+ }
+ }
+ else
+ {
+ //Read the generic_kernel_parameters line
+ fgets(cmd, 80, fp);
+ }
+
+ }
+ else if (strcmp(cmd, "composed_kernel_parameters") == 0)
+ {
+ if (param.kernel_type == COMPOSED)
+ {
+ //Load generic parameters
+ delete generic_kernel_functor;
+
+ param.kernel_composed = new ComposedKernelFunctor;
+ int cr = param.kernel_composed->load_parameters(&fp);
+ model->delete_composed = true;
+ //int cr = param.kernel_generic->load_parameters(&fp);
+ if (cr != 0)
+ {
+ fprintf(stderr, "error while loading composed kernel parameters from the file %s.\n", model_file_name);
+ }
+ }
+ }
+ /*** End OTB modification ***/
+ else
+ {
+ fprintf(stderr,"unknown text in model file: [%s]\n",cmd);
+ free(model->rho);
+ free(model->label);
+ free(model->nSV);
+ free(model);
+ return NULL;
+ }
+ }
+
+ // read sv_coef and SV
+
+ int elements = 0;
+ long pos = ftell(fp);
+
+ max_line_len = 1024;
+ line = Malloc(char,max_line_len);
+ char *p,*endptr,*idx,*val;
+
+ while(readline(fp)!=NULL)
+ {
+ p = strtok(line,":");
+ while(1)
+ {
+ p = strtok(NULL,":");
+ if(p == NULL)
+ break;
+ ++elements;
+ }
+ }
+ elements += model->l;
+
+ fseek(fp,pos,SEEK_SET);
+
+ int m = model->nr_class - 1;
+ int l = model->l;
+ model->sv_coef = Malloc(double *,m);
+ int i;
+ for(i=0;i<m;i++)
+ model->sv_coef[i] = Malloc(double,l);
+ model->SV = Malloc(svm_node*,l);
+ svm_node *x_space = NULL;
+ if(l>0) x_space = Malloc(svm_node,elements);
+
+ int j=0;
+ for(i=0;i<l;i++)
+ {
+ readline(fp);
+ model->SV[i] = &x_space[j];
+
+ p = strtok(line, " \t");
+ model->sv_coef[0][i] = strtod(p,&endptr);
+ for(int k=1;k<m;k++)
+ {
+ p = strtok(NULL, " \t");
+ model->sv_coef[k][i] = strtod(p,&endptr);
+ }
+
+ while(1)
+ {
+ idx = strtok(NULL, ":");
+ val = strtok(NULL, " \t");
+
+ if(val == NULL)
+ break;
+ x_space[j].index = (int) strtol(idx,&endptr,10);
+ x_space[j].value = strtod(val,&endptr);
+
+ ++j;
+ }
+ x_space[j++].index = -1;
+ }
+ free(line);
+
+ if (ferror(fp) != 0 || fclose(fp) != 0)
+ return NULL;
+
+ model->free_sv = 1; // XXX
+ return model;
+}
+
+/*** Begin OTB modification ***/
+svm_model *svm_copy_model( const svm_model *model )
+{
+ const svm_parameter& param = model->param;
+
+ // instanciated the copy
+ svm_model *modelCpy = Malloc(svm_model,1);
+ svm_parameter& paramCpy = modelCpy->param;
+ modelCpy->rho = NULL;
+ modelCpy->probA = NULL;
+ modelCpy->probB = NULL;
+ modelCpy->label = NULL;
+ modelCpy->nSV = NULL;
+ modelCpy->delete_composed = false;
+
+ // SVM type copy
+ paramCpy.svm_type = param.svm_type;
+ // Kernel type copy
+ paramCpy.kernel_type = param.kernel_type;
+ // Param copy
+ paramCpy.degree = param.degree;
+ paramCpy.gamma = param.gamma;
+ paramCpy.coef0 = param.coef0;
+ // Model variable
+ int nr_class = model->nr_class;
+ int l = model->l;
+
+ modelCpy->nr_class = nr_class;
+ modelCpy->l = l;
+ if (model->rho)
+ {
+ int n = model->nr_class * (model->nr_class - 1) / 2;
+ modelCpy->rho = Malloc(double,n);
+ for (int i = 0; i < n; i++)
+ modelCpy->rho[i] = model->rho[i];
+ }
+ if (model->label)
+ {
+ modelCpy->label = Malloc(int,nr_class);
+ for (int i = 0; i < nr_class; i++)
+ modelCpy->label[i] = model->label[i];
+ }
+ if (model->probA)
+ {
+ int n = nr_class * (nr_class - 1) / 2;
+ modelCpy->probA = Malloc(double,n);
+ for (int i = 0; i < n; i++)
+ modelCpy->probA[i] = model->probA[i];
+ }
+ if (model->probB)
+ {
+ int n = nr_class * (nr_class - 1) / 2;
+ modelCpy->probB = Malloc(double,n);
+ for (int i = 0; i < n; i++)
+ modelCpy->probB[i] = model->probB[i];
+ }
+ if (model->nSV)
+ {
+ modelCpy->nSV = Malloc(int,nr_class);
+ for (int i = 0; i < nr_class; i++)
+ modelCpy->nSV[i] = model->nSV[i];
+ }
+
+ // SV copy
+ const double * const *sv_coef = model->sv_coef;
+ const svm_node * const *SV = model->SV;
+
+ modelCpy->SV = Malloc(svm_node*,l);
+ svm_node **SVCpy = modelCpy->SV;
+
+ modelCpy->sv_coef = Malloc(double *,nr_class-1);
+
+ for (int i = 0; i < nr_class - 1; i++)
+ modelCpy->sv_coef[i] = Malloc(double,l);
+
+ // Compute the total number of SV elements.
+ unsigned int elements = 0;
+ for (int p = 0; p < l; p++)
+ {
+ const svm_node *tempNode = SV[p];
+ while (tempNode->index != -1)
+ {
+ tempNode++;
+ elements++;
+ }
+ elements++;// for -1 values
+ }
+
+ if (l > 0)
+ {
+ modelCpy->SV[0] = Malloc(svm_node,elements);
+ memcpy(modelCpy->SV[0], model->SV[0], sizeof(svm_node*) * elements);
+ }
+ svm_node *x_space = modelCpy->SV[0];
+
+ int j = 0;
+ for (int i = 0; i < l; i++)
+ {
+ // sv_coef
+ for (int k = 0; k < nr_class - 1; k++)
+ modelCpy->sv_coef[k][i] = sv_coef[k][i];
+
+ // SV
+ modelCpy->SV[i] = &x_space[j];
+ const svm_node *p = SV[i];
+ svm_node *pCpy = SVCpy[i];
+ while (p->index != -1)
+ {
+ pCpy->index = p->index;
+ pCpy->value = p->value;
+ p++;
+ pCpy++;
+ j++;
+ }
+ pCpy->index = -1;
+ j++;
+ }
+
+ // Generic kernel copy
+ if (param.kernel_type == GENERIC)
+ {
+ paramCpy.kernel_generic = Malloc(GenericKernelFunctorBase, sizeof(*(param.kernel_generic)));
+ memcpy(paramCpy.kernel_generic, param.kernel_generic, sizeof(*(param.kernel_generic)));
+ }
+ // Composrd kernel copy
+ if (param.kernel_type == COMPOSED)
+ {
+ paramCpy.kernel_composed = Malloc(ComposedKernelFunctor, 1);
+ *(paramCpy.kernel_composed) = *(param.kernel_composed);
+ }
+
+ return modelCpy;
+}
+/*** End OTB modification ***/
+
+
+
+void svm_free_model_content(svm_model* model_ptr)
+{
+ if(model_ptr->free_sv && model_ptr->l > 0)
+ free((void *)(model_ptr->SV[0]));
+ for(int i=0;i<model_ptr->nr_class-1;i++)
+ free(model_ptr->sv_coef[i]);
+ free(model_ptr->SV);
+ free(model_ptr->sv_coef);
+ free(model_ptr->rho);
+ free(model_ptr->label);
+ free(model_ptr->probA);
+ free(model_ptr->probB);
+ free(model_ptr->nSV);
+}
+
+void svm_free_and_destroy_model(svm_model** model_ptr_ptr)
+{
+ svm_model* model_ptr = *model_ptr_ptr;
+ if(model_ptr != NULL)
+ {
+ svm_free_model_content(model_ptr);
+ free(model_ptr);
+ }
+}
+
+void svm_destroy_model(svm_model* model_ptr)
+{
+ fprintf(stderr,"warning: svm_destroy_model is deprecated and should not be used. Please use svm_free_and_destroy_model(svm_model **model_ptr_ptr)\n");
+ svm_free_and_destroy_model(&model_ptr);
+}
+
+void svm_destroy_param(svm_parameter* param)
+{
+ free(param->weight_label);
+ free(param->weight);
+}
+
+const char *svm_check_parameter(const svm_problem *prob, const svm_parameter *param)
+{
+ // svm_type
+
+ int svm_type = param->svm_type;
+ if(svm_type != C_SVC &&
+ svm_type != NU_SVC &&
+ svm_type != ONE_CLASS &&
+ svm_type != EPSILON_SVR &&
+ svm_type != NU_SVR)
+ return "unknown svm type";
+
+ // kernel_type, degree
+
+ int kernel_type = param->kernel_type;
+ if(kernel_type != LINEAR &&
+ kernel_type != POLY &&
+ kernel_type != RBF &&
+ kernel_type != SIGMOID &&
+ kernel_type != PRECOMPUTED
+ /*** Begin OTB modification ***/
+ && kernel_type != GENERIC
+ && kernel_type != COMPOSED
+ /*** End OTB modification ***/
+ )
+ return "unknown kernel type";
+
+ if(param->gamma < 0)
+ return "gamma < 0";
+
+ if(param->degree < 0)
+ return "degree of polynomial kernel < 0";
+
+ // cache_size,eps,C,nu,p,shrinking
+
+ if(param->cache_size <= 0)
+ return "cache_size <= 0";
+
+ if(param->eps <= 0)
+ return "eps <= 0";
+
+ if(svm_type == C_SVC ||
+ svm_type == EPSILON_SVR ||
+ svm_type == NU_SVR)
+ if(param->C <= 0)
+ return "C <= 0";
+
+ if(svm_type == NU_SVC ||
+ svm_type == ONE_CLASS ||
+ svm_type == NU_SVR)
+ if(param->nu <= 0 || param->nu > 1)
+ return "nu <= 0 or nu > 1";
+
+ if(svm_type == EPSILON_SVR)
+ if(param->p < 0)
+ return "p < 0";
+
+ if(param->shrinking != 0 &&
+ param->shrinking != 1)
+ return "shrinking != 0 and shrinking != 1";
+
+ if(param->probability != 0 &&
+ param->probability != 1)
+ return "probability != 0 and probability != 1";
+
+ if(param->probability == 1 &&
+ svm_type == ONE_CLASS)
+ return "one-class SVM probability output not supported yet";
+
+
+ // check whether nu-svc is feasible
+
+ if(svm_type == NU_SVC)
+ {
+ int l = prob->l;
+ int max_nr_class = 16;
+ int nr_class = 0;
+ int *label = Malloc(int,max_nr_class);
+ int *count = Malloc(int,max_nr_class);
+
+ int i;
+ for(i=0;i<l;i++)
+ {
+ int this_label = (int)prob->y[i];
+ int j;
+ for(j=0;j<nr_class;j++)
+ if(this_label == label[j])
+ {
+ ++count[j];
+ break;
+ }
+ if(j == nr_class)
+ {
+ if(nr_class == max_nr_class)
+ {
+ max_nr_class *= 2;
+ label = (int *)realloc(label,max_nr_class*sizeof(int));
+ count = (int *)realloc(count,max_nr_class*sizeof(int));
+ }
+ label[nr_class] = this_label;
+ count[nr_class] = 1;
+ ++nr_class;
+ }
+ }
+
+ for(i=0;i<nr_class;i++)
+ {
+ int n1 = count[i];
+ for(int j=i+1;j<nr_class;j++)
+ {
+ int n2 = count[j];
+ if(param->nu*(n1+n2)/2 > min(n1,n2))
+ {
+ free(label);
+ free(count);
+ return "specified nu is infeasible";
+ }
+ }
+ }
+ free(label);
+ free(count);
+ }
+
+ return NULL;
+}
+
+int svm_check_probability_model(const svm_model *model)
+{
+ return ((model->param.svm_type == C_SVC || model->param.svm_type == NU_SVC) &&
+ model->probA!=NULL && model->probB!=NULL) ||
+ ((model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) &&
+ model->probA!=NULL);
+}
+
+void svm_set_print_string_function(void (*print_func)(const char *))
+{
+ if(print_func == NULL)
+ svm_print_string = &print_string_stdout;
+ else
+ svm_print_string = print_func;
+}
+
+/*** Begin OTB modification ***/
+
+GenericKernelFunctorBase::GenericKernelFunctorBase(const GenericKernelFunctorBase& copy)
+{
+ *this = copy;
+}
+
+GenericKernelFunctorBase&
+GenericKernelFunctorBase::operator=(const GenericKernelFunctorBase& copy)
+{
+ this->m_MapParameters = copy.m_MapParameters;
+ this->m_Name = copy.m_Name;
+ return *this;
+}
+
+int
+GenericKernelFunctorBase::
+load_parameters(FILE ** pfile)
+{
+ int NbParams(0);
+ char keyword[81];
+ char value[81];
+ // Read functor name
+ fscanf(*pfile,"%80s",keyword);
+ m_Name = std::string(keyword);
+ // Read number of parameters
+ fscanf(*pfile,"%d",&NbParams);
+ // if( NbParams == 0 ) return -1;
+ for ( int cpt=0 ; cpt < NbParams ; cpt++)
+ {
+ fscanf(*pfile,"%80s",keyword);
+ fscanf(*pfile,"%80s",value);
+ m_MapParameters[std::string(keyword)] = std::string(value);
+ }
+ return 0;
+}
+int
+GenericKernelFunctorBase::
+save_parameters(FILE ** pfile, const char * generic_kernel_parameters_keyword)const
+{
+ MapConstIterator iter=m_MapParameters.begin();
+
+ std::string line(generic_kernel_parameters_keyword);
+ std::string strNbParams;
+ ::otb::StringStream flux;
+ flux << m_MapParameters.size();
+ flux >> strNbParams;
+ line = line + " " + m_Name + " " + strNbParams;
+
+ // DON'T USE ITER because for a COPY ONLY, whereas map.size() IS GOOD, map.end()
+ // doesn't give the good answer -> SEGFAULT
+ // while( iter != m_MapParameters.end() )
+ // {
+ // line = line + " " + iter->first + " " + iter->second;
+ // ++iter;
+ // }
+ for ( unsigned int i = 0; i<m_MapParameters.size(); i++)
+ {
+ line = line + " " + iter->first + " " + iter->second;
+ ++iter;
+ }
+ fprintf(*pfile,"%s\n", line.c_str());
+ return 0;
+}
+
+void
+GenericKernelFunctorBase::
+print_parameters(void)const
+{
+ MapConstIterator iter=this->GetMapParameters().begin();
+
+ std::cout << "Print generic kernel parameters: "<<this->GetMapParameters().size()<<std::endl;
+ while( iter != this->GetMapParameters().end() )
+ {
+ std::cout << " "<<iter->first <<" "<<iter->second<<std::endl;
+ ++iter;
+ }
+}
+
+double
+GenericKernelFunctorBase::
+dot(const svm_node *px, const svm_node *py)const
+{
+ double sum = 0.;
+
+ while(px->index != -1 && py->index != -1)
+ {
+ if(px->index == py->index)
+ {
+ sum += px->value * py->value;
+ ++px;
+ ++py;
+ }
+ else
+ {
+ if(px->index > py->index)
+ ++py;
+ else
+ ++px;
+ }
+ }
+ return sum;
+}
+
+svm_node *
+GenericKernelFunctorBase::
+sub(const svm_node *px, const svm_node *py) const
+ /* compute the difference a-b of two sparse vectors */
+ /* Note: SVECTOR lists are not followed, but only the first
+ SVECTOR is used */
+{
+ long veclength = 1;
+
+ const svm_node * pxbis = px;
+ const svm_node * pybis = py;
+
+ while (px->index != -1 && py->index != -1)
+ {
+ ++veclength;
+ if(px->index == py->index)
+ {
+ ++px;
+ ++py;
+ }
+ else
+ {
+ if (px->index < py->index)
+ {
+ ++px;
+ }
+ else
+ {
+ ++py;
+ }
+ }
+ }
+
+ while (py->index != -1)
+ {
+ ++veclength;
+ ++py;
+ }
+ while (px->index != -1)
+ {
+ ++veclength;
+ ++px;
+ }
+
+ svm_node *vec;
+ vec = new svm_node[veclength];
+ unsigned long int vecIt = 0;
+
+ px = pxbis;
+ py = pybis;
+
+ while (px->index != -1 && py->index != -1)
+ {
+ if(px->index == py->index)
+ {
+ (vec[vecIt])=(*px);
+ vec[vecIt].value-=py->value;
+
+ if(vec[vecIt].value != 0)
+ {
+ ++vecIt;
+ }
+ ++px;
+ ++py;
+
+ }
+ else
+ {
+ if(px->index > py->index)
+ {
+ (vec[vecIt])=(*py);
+ vec[vecIt].value*=(-1);
+ ++vecIt;
+ ++py;
+ }
+ else
+ {
+ (vec[vecIt])=(*px);
+ ++vecIt;
+ ++px;
+ }
+ }
+ }
+
+ while (py->index != -1)
+ {
+ (vec[vecIt])=(*py);
+ vec[vecIt].value*=(-1);
+ ++vecIt;
+ ++py;
+ }
+ while (px->index != -1)
+ {
+ (vec[vecIt])=(*px);
+ ++vecIt;
+ ++px;
+ }
+ return(vec);
+}
+
+
+svm_node *
+GenericKernelFunctorBase::
+add(const svm_node *px, const svm_node *py) const
+ /* compute the sum a+b of two sparse vectors */
+ /* Note: SVECTOR lists are not followed, but only the first
+ SVECTOR is used */
+{
+ const svm_node * pxbis = px;
+ const svm_node * pybis = py;
+ long veclength = 1;
+
+ while (px->index != -1 && py->index != -1)
+ {
+ ++veclength;
+ if(px->index == py->index)
+ {
+ ++px;
+ ++py;
+ }
+ else
+ {
+ if (px->index < py->index)
+ {
+ ++px;
+ }
+ else
+ {
+ ++py;
+ }
+ }
+ }
+
+ while (py->index != -1)
+ {
+ ++veclength;
+ ++py;
+ }
+ while (px->index != -1)
+ {
+ ++veclength;
+ ++px;
+ }
+
+ svm_node *vec;
+ vec = new svm_node[veclength];
+ unsigned long int vecIt = 0;
+
+ px = pxbis;
+ py = pybis;
+
+ while (px->index != -1 && py->index != -1)
+ {
+ if(px->index == py->index)
+ {
+ (vec[vecIt])=(*px);
+ vec[vecIt].value+=py->value;
+ if(vec[vecIt].value != 0)
+ {
+ ++vecIt;
+ }
+ ++px;
+ ++py;
+ }
+ else
+ {
+ if (px->index < py->index)
+ {
+ (vec[vecIt])=(*px);
+ ++vecIt;
+ ++px;
+ }
+ else
+ {
+ (vec[vecIt])=(*py);
+ ++vecIt;
+ ++py;
+
+ }
+ }
+ }
+ while (py->index != -1)
+ {
+ (vec[vecIt])=(*py);
+ ++vecIt;
+ ++py;
+ }
+ while (px->index != -1)
+ {
+ (vec[vecIt])=(*px);
+ ++vecIt;
+ ++px;
+ }
+
+ return(vec);
+}
+
+// ****************************************************************************************
+// ************************ ComposedKernelFunctor methods ********************/
+// ****************************************************************************************
+
+ComposedKernelFunctor::ComposedKernelFunctor(const ComposedKernelFunctor& copy)
+{
+ *this = copy;
+}
+
+ComposedKernelFunctor&
+ComposedKernelFunctor::operator=(const ComposedKernelFunctor& copy)
+{
+ // Call Superclass::operator=
+ static_cast<Superclass&>(*this) = static_cast<const Superclass&>(copy);
+
+ // Copy Self attributes
+ this->m_KernelFunctorList = copy.m_KernelFunctorList;
+ this->m_HaveToBeDeletedList = copy.m_HaveToBeDeletedList;
+ this->m_PonderationList = copy.m_PonderationList;
+ return *this;
+}
+
+void
+ComposedKernelFunctor
+::print_parameters(void)const
+{
+ MapConstIterator iter = this->GetMapParameters().begin();
+ std::cout << "Print composed kernel parameters: "<<this->GetName()<<", "<<this->GetMapParameters().size()<<std::endl;
+ while( iter != this->GetMapParameters().end() )
+ {
+ std::cout << " "<<iter->first <<" "<<iter->second<<std::endl;
+ ++iter;
+ }
+ std::cout<<std::endl;
+ std::cout<<"Composition kernels:"<<std::endl;
+ if (m_KernelFunctorList.size() != 0 && m_PonderationList.size() != 0 && m_KernelFunctorList.size() == m_PonderationList.size())
+ {
+ for (unsigned int i = 0; i<m_KernelFunctorList.size(); i++)
+ {
+ std::cout<<m_KernelFunctorList[i]->GetName()<<":"<<std::endl;
+ std::cout<<"Associated ponderation:"<<m_PonderationList[i]<<std::endl;
+ m_KernelFunctorList[i]->print_parameters();
+ std::cout<<std::endl;
+ }
+ }
+ else
+ {
+ itkGenericExceptionMacro(<<"ComposedKernelFunctor::print_param() : lists dimensions mismatch");
+ }
+}
+
+
+int
+ComposedKernelFunctor::
+load_parameters(FILE ** pfile)
+{
+ int NbParams(0);
+ char keyword[81];
+ char value[81];
+
+ // Read functor name
+ fscanf(*pfile,"%80s",keyword);
+ this->SetName(std::string(keyword));
+ // Read number of parameters
+ fscanf(*pfile,"%d",&NbParams);
+
+ for ( int cpt=0 ; cpt < NbParams ; cpt++)
+ {
+ fscanf(*pfile,"%80s",keyword);
+ fscanf(*pfile,"%80s",value);
+ this->SetValue<std::string>(keyword, value);
+ }
+
+ char tempChar[100];
+ fscanf(*pfile, "%80s", tempChar);
+ while( strcmp(tempChar,"Ponderation")==0 || strcmp(tempChar,"list:")==0 )
+ {
+ fscanf(*pfile, "%80s", tempChar);
+ }
+
+ unsigned int i = 0;
+ while( strcmp(tempChar,"Kernels")!=0 )
+ {
+ m_PonderationList.push_back(::atof(tempChar));
+ fscanf(*pfile, "%80s", tempChar);
+ i++;
+ }
+ while( strcmp(tempChar,"Kernels")==0 || strcmp(tempChar,"list:")==0 || strcmp(tempChar,"Number")==0 || strcmp(tempChar,"of")==0 || strcmp(tempChar,"Kernels:")==0)
+ {
+ fscanf(*pfile, "%80s", tempChar);
+ }
+ int NbOfKernels = ::atoi(tempChar);
+
+ for(unsigned int j=0; j<static_cast<unsigned int>(NbOfKernels); j++)
+ {
+ fscanf(*pfile, "%80s", tempChar);
+ GenericKernelFunctorBase * gen;
+ gen = new GenericKernelFunctorBase;
+ gen->load_parameters(pfile);
+ m_KernelFunctorList.push_back(gen);
+ // Add the pointer to the "Have to Deleted" pointer list
+ m_HaveToBeDeletedList.push_back(gen);
+ }
+
+ return 0;
+}
+
+
+int
+ComposedKernelFunctor::
+save_parameters(FILE ** pfile, const char * composed_kernel_parameters_keyword)const
+{
+ MapConstIterator iter = this->GetMapParameters().begin();
+ std::string line(composed_kernel_parameters_keyword);
+ std::string strNbParams;
+ ::otb::StringStream flux;
+ flux << this->GetMapParameters().size();
+ flux >> strNbParams;
+ line = line + " " + this->GetName() + " " + strNbParams;
+ while( iter != this->GetMapParameters().end() )
+ {
+ line = line + " " + iter->first + " " + iter->second;
+ ++iter;
+ }
+ line = line + "\n" + "Ponderation list:\n";
+
+ for (unsigned int i = 0; i<m_PonderationList.size(); i++)
+ {
+ std::string ponde;
+ ::otb::StringStream flux;
+ flux << m_PonderationList[i];
+ flux >> ponde;
+ line = line + " " + ponde;
+ }
+ line = line + "\n" + "Kernels list:\nNumber of Kernels: ";
+ std::string nbOfKernels;
+ ::otb::StringStream flux2;
+ flux2 << m_KernelFunctorList.size();
+ flux2 >> nbOfKernels;
+
+ line = line + nbOfKernels + "\n";
+ fprintf(*pfile,"%s", line.c_str());
+ for (unsigned int i = 0; i<m_KernelFunctorList.size(); i++)
+ {
+ m_KernelFunctorList[i]->save_parameters(pfile, "generic_kernel_parameters");
+ }
+
+ return 0;
+}
+/*** End OTB modification ***/
+
+
diff --git a/Utilities/otbsvm/svm.cxx b/Utilities/otbsvm/svm.cxx
deleted file mode 100644
index dd5ae3d2d5..0000000000
--- a/Utilities/otbsvm/svm.cxx
+++ /dev/null
@@ -1,3788 +0,0 @@
-//OTB's modifications
-#include <iostream>
-
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <ctype.h>
-#include <float.h>
-#include <string.h>
-#include <stdarg.h>
-#include "svm.h"
-typedef float Qfloat;
-//#include <double.h>
-//typedef double Qfloat;
-
-
-typedef signed char schar;
-#ifndef min
-template <class T> inline T min(T x,T y) { return (x<y)?x:y; }
-#endif
-#ifndef max
-template <class T> inline T max(T x,T y) { return (x>y)?x:y; }
-#endif
-template <class T> inline void swap(T& x, T& y) { T t=x; x=y; y=t; }
-template <class S, class T> inline void clone(T*& dst, S* src, int n)
-{
- dst = new T[n];
- memcpy((void *)dst,(void *)src,sizeof(T)*n);
-}
-inline double powi(double base, int times)
-{
- double tmp = base, ret = 1.0;
-
- for(int t=times; t>0; t/=2)
- {
- if(t%2==1) ret*=tmp;
- tmp = tmp * tmp;
- }
- return ret;
-}
-#define INF HUGE_VAL
-#define TAU 1e-12
-#define Malloc(type,n) new type[n]
-#if 1
-void info(char *fmt,...)
-{
- va_list ap;
- va_start(ap,fmt);
- vprintf(fmt,ap);
- va_end(ap);
-}
-void info_flush()
-{
- fflush(stdout);
-}
-#else
-void info(char *fmt,...) {}
-void info_flush() {}
-#endif
-
-// ****************************************************************************************
-// Kernel Cache
-// ****************************************************************************************
-// l is the number of total data items
-// size is the cache size limit in bytes
-//
-class Cache
-{
-public:
- Cache(int l,long int size);
- ~Cache();
-
- // request data [0,len)
- // return some position p where [p,len) need to be filled
- // (p >= len if nothing needs to be filled)
- int get_data(const int index, Qfloat **data, int len);
- void swap_index(int i, int j); // future_option
-private:
- int l;
- long int size;
- struct head_t
- {
- head_t *prev, *next; // a cicular list
- Qfloat *data;
- int len; // data[0,len) is cached in this entry
- };
-
- head_t *head;
- head_t lru_head;
- void lru_delete(head_t *h);
- void lru_insert(head_t *h);
-};
-
-Cache::Cache(int l_,long int size_):l(l_),size(size_)
-{
- head = (head_t *)calloc(l,sizeof(head_t)); // initialized to 0
- size /= sizeof(Qfloat);
- size -= l * sizeof(head_t) / sizeof(Qfloat);
- size = max(size, (long int) 2*l); // cache must be large enough for two columns
- lru_head.next = lru_head.prev = &lru_head;
-}
-
-Cache::~Cache()
-{
- for(head_t *h = lru_head.next; h != &lru_head; h=h->next)
- delete [](h->data);
- delete [](head);
-}
-
-void Cache::lru_delete(head_t *h)
-{
- // delete from current location
- h->prev->next = h->next;
- h->next->prev = h->prev;
-}
-
-void Cache::lru_insert(head_t *h)
-{
- // insert to last position
- h->next = &lru_head;
- h->prev = lru_head.prev;
- h->prev->next = h;
- h->next->prev = h;
-}
-
-int Cache::get_data(const int index, Qfloat **data, int len)
-{
- head_t *h = &head[index];
- if(h->len) lru_delete(h);
- int more = len - h->len;
-
- if(more > 0)
- {
- // free old space
- while(size < more)
- {
- head_t *old = lru_head.next;
- lru_delete(old);
- delete [](old->data);
- size += old->len;
- old->data = 0;
- old->len = 0;
- }
-
- // allocate new space
- h->data = (Qfloat *)realloc(h->data,sizeof(Qfloat)*len);
- size -= more;
- swap(h->len,len);
- }
-
- lru_insert(h);
- *data = h->data;
- return len;
-}
-
-void Cache::swap_index(int i, int j)
-{
- if(i==j) return;
-
- if(head[i].len) lru_delete(&head[i]);
- if(head[j].len) lru_delete(&head[j]);
- swap(head[i].data,head[j].data);
- swap(head[i].len,head[j].len);
- if(head[i].len) lru_insert(&head[i]);
- if(head[j].len) lru_insert(&head[j]);
-
- if(i>j) swap(i,j);
- for(head_t *h = lru_head.next; h!=&lru_head; h=h->next)
- {
- if(h->len > i)
- {
- if(h->len > j)
- swap(h->data[i],h->data[j]);
- else
- {
- // give up
- lru_delete(h);
- delete [](h->data);
- size += h->len;
- h->data = 0;
- h->len = 0;
- }
- }
- }
-}
-
-
-// ****************************************************************************************
-// Kernel evaluation
-// ****************************************************************************************
-// the static method k_function is for doing single kernel evaluation
-// the constructor of Kernel prepares to calculate the l*l kernel matrix
-// the member function get_Q is for getting one column from the Q Matrix
-//
-class QMatrix {
-public:
- virtual Qfloat *get_Q(int column, int len) const = 0;
- virtual Qfloat *get_QD() const = 0;
- virtual void swap_index(int i, int j) const = 0;
- virtual ~QMatrix() {}
-};
-
-class Kernel: public QMatrix {
-public:
- Kernel(int l, svm_node * const * x, const svm_parameter& param);
- virtual ~Kernel();
-
- static double k_function(const svm_node *x, const svm_node *y,
- const svm_parameter& param);
- virtual Qfloat *get_Q(int column, int len) const = 0;
- virtual Qfloat *get_QD() const = 0;
- virtual void swap_index(int i, int j) const // no so const...
- {
- swap(x[i],x[j]);
- if(x_square) swap(x_square[i],x_square[j]);
- }
-
-protected:
-
- //OTB's modifications
- double (Kernel::*kernel_function)(int i, int j, const svm_parameter& param) const;
- const svm_parameter& m_param;
-
-private:
- const svm_node **x;
- double *x_square;
-
- // svm_parameter
- const int kernel_type;
- const int degree;
- const double gamma;
- const double coef0;
-
- static double dot(const svm_node *px, const svm_node *py);
- double kernel_linear(int i, int j, const svm_parameter& param) const
- {
- return dot(x[i],x[j]);
- }
- double kernel_poly(int i, int j, const svm_parameter& param) const
- {
- return powi(gamma*dot(x[i],x[j])+coef0,degree);
- }
- double kernel_rbf(int i, int j, const svm_parameter& param) const
- {
- return exp(-gamma*(x_square[i]+x_square[j]-2*dot(x[i],x[j])));
- }
- double kernel_sigmoid(int i, int j, const svm_parameter& param) const
- {
- return tanh(gamma*dot(x[i],x[j])+coef0);
- }
- double kernel_precomputed(int i, int j, const svm_parameter& param) const
- {
- return x[i][(int)(x[j][0].value)].value;
- }
- //OTB's modifications
- double kernel_generic(int i, int j, const svm_parameter& param) const
- {
- if( param.kernel_generic == NULL )
- {
- itkGenericExceptionMacro( << "Generic Kernel is not initialiszed !");
- }
- return ((*param.kernel_generic)(x[i],x[j],param));
- }
- // Composed kernel
- double kernel_composed(int i, int j, const svm_parameter& param) const
- {
- if( param.kernel_composed == NULL )
- {
- itkGenericExceptionMacro( << "Generic Kernel is not initialiszed !");
- }
- return ((*param.kernel_composed)(x[i],x[j],param));
- }
-};
-
-Kernel::Kernel(int l, svm_node * const * x_, const svm_parameter& param)
-//OTB's modifications
-:m_param(param),kernel_type(param.kernel_type), degree(param.degree),
- gamma(param.gamma), coef0(param.coef0)
-{
- switch(kernel_type)
- {
- case LINEAR:
- kernel_function = &Kernel::kernel_linear;
- break;
- case POLY:
- kernel_function = &Kernel::kernel_poly;
- break;
- case RBF:
- kernel_function = &Kernel::kernel_rbf;
- break;
- case SIGMOID:
- kernel_function = &Kernel::kernel_sigmoid;
- break;
- case PRECOMPUTED:
- kernel_function = &Kernel::kernel_precomputed;
- break;
- //OTB's modifications
- case GENERIC:
- kernel_function = &Kernel::kernel_generic;
- break;
- case COMPOSED:
- kernel_function = &Kernel::kernel_composed;
- break;
-
- }
-
- clone(x,x_,l);
-
- if(kernel_type == RBF)
- {
- x_square = new double[l];
- for(int i=0;i<l;i++)
- x_square[i] = dot(x[i],x[i]);
- }
- else
- x_square = 0;
-}
-
-Kernel::~Kernel()
-{
- delete[] x;
- delete[] x_square;
-}
-
-double Kernel::dot(const svm_node *px, const svm_node *py)
-{
- double sum = 0.;
- while(px->index != -1 && py->index != -1)
- {
- if(px->index == py->index)
- {
- sum += px->value * py->value;
- ++px;
- ++py;
- }
- else
- {
- if(px->index > py->index)
- ++py;
- else
- ++px;
- }
- }
- return sum;
-}
-
-double Kernel::k_function(const svm_node *x, const svm_node *y,
- const svm_parameter& param)
-{
- switch(param.kernel_type)
- {
- case LINEAR:
- return dot(x,y);
- case POLY:
- return powi(param.gamma*dot(x,y)+param.coef0,param.degree);
- case RBF:
- {
- double sum = 0;
- while(x->index != -1 && y->index !=-1)
- {
- if(x->index == y->index)
- {
- double d = x->value - y->value;
- sum += d*d;
- ++x;
- ++y;
- }
- else
- {
- if(x->index > y->index)
- {
- sum += y->value * y->value;
- ++y;
- }
- else
- {
- sum += x->value * x->value;
- ++x;
- }
- }
- }
-
- while(x->index != -1)
- {
- sum += x->value * x->value;
- ++x;
- }
-
- while(y->index != -1)
- {
- sum += y->value * y->value;
- ++y;
- }
-
- return exp(-param.gamma*sum);
- }
- case SIGMOID:
- return tanh(param.gamma*dot(x,y)+param.coef0);
- case PRECOMPUTED: //x: test (validation), y: SV
- return x[(int)(y->value)].value;
- //OTB's modifications
- case GENERIC:
- return ((*param.kernel_generic)(x,y,param));
- // Composed kernel
- case COMPOSED:
- return ((*param.kernel_composed)(x,y,param));
- default:
- return 0; /* Unreachable */
- }
-}
-
-
-// ****************************************************************************************
-// SOLVER
-// ****************************************************************************************
-// An SMO algorithm in Fan et al., JMLR 6(2005), p. 1889--1918
-// Solves:
-//
-// min 0.5(\alpha^T Q \alpha) + p^T \alpha
-//
-// y^T \alpha = \delta
-// y_i = +1 or -1
-// 0 <= alpha_i <= Cp for y_i = 1
-// 0 <= alpha_i <= Cn for y_i = -1
-//
-// Given:
-//
-// Q, p, y, Cp, Cn, and an initial feasible point \alpha
-// l is the size of vectors and matrices
-// eps is the stopping tolerance
-//
-// solution will be put in \alpha, objective value will be put in obj
-//
-class Solver {
-public:
- Solver() {};
- virtual ~Solver() {};
-
- struct SolutionInfo {
- double obj;
- double rho;
- double upper_bound_p;
- double upper_bound_n;
- double r; // for Solver_NU
- };
-
- void Solve(int l, const QMatrix& Q, const double *p_, const schar *y_,
- double *alpha_, double Cp, double Cn, double eps,
- SolutionInfo* si, int shrinking);
-protected:
- int active_size;
- schar *y;
- double *G; // gradient of objective function
- enum { LOWER_BOUND, UPPER_BOUND, FREE };
- char *alpha_status; // LOWER_BOUND, UPPER_BOUND, FREE
- double *alpha;
- const QMatrix *Q;
- const Qfloat *QD;
- double eps;
- double Cp,Cn;
- double *p;
- int *active_set;
- double *G_bar; // gradient, if we treat free variables as 0
- int l;
- bool unshrinked; // XXX
-
- double get_C(int i)
- {
- return (y[i] > 0)? Cp : Cn;
- }
- void update_alpha_status(int i)
- {
- if(alpha[i] >= get_C(i))
- alpha_status[i] = UPPER_BOUND;
- else if(alpha[i] <= 0)
- alpha_status[i] = LOWER_BOUND;
- else alpha_status[i] = FREE;
- }
- bool is_upper_bound(int i) { return alpha_status[i] == UPPER_BOUND; }
- bool is_lower_bound(int i) { return alpha_status[i] == LOWER_BOUND; }
- bool is_free(int i) { return alpha_status[i] == FREE; }
- void swap_index(int i, int j);
- void reconstruct_gradient();
- virtual int select_working_set(int &i, int &j);
- virtual double calculate_rho();
- virtual void do_shrinking();
-private:
- bool be_shrunken(int i, double Gmax1, double Gmax2);
-};
-
-void Solver::swap_index(int i, int j)
-{
- Q->swap_index(i,j);
- swap(y[i],y[j]);
- swap(G[i],G[j]);
- swap(alpha_status[i],alpha_status[j]);
- swap(alpha[i],alpha[j]);
- swap(p[i],p[j]);
- swap(active_set[i],active_set[j]);
- swap(G_bar[i],G_bar[j]);
-}
-
-void Solver::reconstruct_gradient()
-{
- // reconstruct inactive elements of G from G_bar and free variables
-
- if(active_size == l) return;
-
- int i;
- for(i=active_size;i<l;i++)
- G[i] = G_bar[i] + p[i];
-
- for(i=0;i<active_size;i++)
- if(is_free(i))
- {
- const Qfloat *Q_i = Q->get_Q(i,l);
- double alpha_i = alpha[i];
- for(int j=active_size;j<l;j++)
- G[j] += alpha_i * Q_i[j];
- }
-}
-
-void Solver::Solve(int l, const QMatrix& Q, const double *p_, const schar *y_,
- double *alpha_, double Cp, double Cn, double eps,
- SolutionInfo* si, int shrinking)
-{
- this->l = l;
- this->Q = &Q;
- QD=Q.get_QD();
- clone(p, p_,l);
- clone(y, y_,l);
- clone(alpha,alpha_,l);
- this->Cp = Cp;
- this->Cn = Cn;
- this->eps = eps;
- unshrinked = false;
-
- // initialize alpha_status
- {
- alpha_status = new char[l];
- for(int i=0;i<l;i++)
- update_alpha_status(i);
- }
-
- // initialize active set (for shrinking)
- {
- active_set = new int[l];
- for(int i=0;i<l;i++)
- active_set[i] = i;
- active_size = l;
- }
-
- // initialize gradient
- {
- G = new double[l];
- G_bar = new double[l];
- int i;
- for(i=0;i<l;i++)
- {
- G[i] = p[i];
- G_bar[i] = 0;
- }
- for(i=0;i<l;i++)
- if(!is_lower_bound(i))
- {
- const Qfloat *Q_i = Q.get_Q(i,l);
- double alpha_i = alpha[i];
- int j;
- for(j=0;j<l;j++)
- G[j] += alpha_i*Q_i[j];
- if(is_upper_bound(i))
- for(j=0;j<l;j++)
- G_bar[j] += get_C(i) * Q_i[j];
- }
- }
-
- // optimization step
-
- int iter = 0;
- int counter = min(l,1000)+1;
-
- while(1)
- {
- // show progress and do shrinking
-
- if(--counter == 0)
- {
- counter = min(l,1000);
- if(shrinking) do_shrinking();
- info("."); info_flush();
- }
-
- int i,j;
- if(select_working_set(i,j)!=0)
- {
- // reconstruct the whole gradient
- reconstruct_gradient();
- // reset active set size and check
- active_size = l;
- info("*"); info_flush();
- if(select_working_set(i,j)!=0)
- break;
- else
- counter = 1; // do shrinking next iteration
- }
-
- ++iter;
-
- // update alpha[i] and alpha[j], handle bounds carefully
-
- const Qfloat *Q_i = Q.get_Q(i,active_size);
- const Qfloat *Q_j = Q.get_Q(j,active_size);
-
- double C_i = get_C(i);
- double C_j = get_C(j);
-
- double old_alpha_i = alpha[i];
- double old_alpha_j = alpha[j];
-
- if(y[i]!=y[j])
- {
- double quad_coef = Q_i[i]+Q_j[j]+2*Q_i[j];
- if (quad_coef <= 0)
- quad_coef = TAU;
- double delta = (-G[i]-G[j])/quad_coef;
- double diff = alpha[i] - alpha[j];
- alpha[i] += delta;
- alpha[j] += delta;
-
- if(diff > 0)
- {
- if(alpha[j] < 0)
- {
- alpha[j] = 0;
- alpha[i] = diff;
- }
- }
- else
- {
- if(alpha[i] < 0)
- {
- alpha[i] = 0;
- alpha[j] = -diff;
- }
- }
- if(diff > C_i - C_j)
- {
- if(alpha[i] > C_i)
- {
- alpha[i] = C_i;
- alpha[j] = C_i - diff;
- }
- }
- else
- {
- if(alpha[j] > C_j)
- {
- alpha[j] = C_j;
- alpha[i] = C_j + diff;
- }
- }
- }
- else
- {
- double quad_coef = Q_i[i]+Q_j[j]-2*Q_i[j];
- if (quad_coef <= 0)
- quad_coef = TAU;
- double delta = (G[i]-G[j])/quad_coef;
- double sum = alpha[i] + alpha[j];
- alpha[i] -= delta;
- alpha[j] += delta;
-
- if(sum > C_i)
- {
- if(alpha[i] > C_i)
- {
- alpha[i] = C_i;
- alpha[j] = sum - C_i;
- }
- }
- else
- {
- if(alpha[j] < 0)
- {
- alpha[j] = 0;
- alpha[i] = sum;
- }
- }
- if(sum > C_j)
- {
- if(alpha[j] > C_j)
- {
- alpha[j] = C_j;
- alpha[i] = sum - C_j;
- }
- }
- else
- {
- if(alpha[i] < 0)
- {
- alpha[i] = 0;
- alpha[j] = sum;
- }
- }
- }
-
- // update G
-
- double delta_alpha_i = alpha[i] - old_alpha_i;
- double delta_alpha_j = alpha[j] - old_alpha_j;
-
- for(int k=0;k<active_size;k++)
- {
- G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j;
- }
-
- // update alpha_status and G_bar
-
- {
- bool ui = is_upper_bound(i);
- bool uj = is_upper_bound(j);
- update_alpha_status(i);
- update_alpha_status(j);
- int k;
- if(ui != is_upper_bound(i))
- {
- Q_i = Q.get_Q(i,l);
- if(ui)
- for(k=0;k<l;k++)
- G_bar[k] -= C_i * Q_i[k];
- else
- for(k=0;k<l;k++)
- G_bar[k] += C_i * Q_i[k];
- }
-
- if(uj != is_upper_bound(j))
- {
- Q_j = Q.get_Q(j,l);
- if(uj)
- for(k=0;k<l;k++)
- G_bar[k] -= C_j * Q_j[k];
- else
- for(k=0;k<l;k++)
- G_bar[k] += C_j * Q_j[k];
- }
- }
- }
-
- // calculate rho
-
- si->rho = calculate_rho();
-
- // calculate objective value
- {
- double v = 0;
- int i;
- for(i=0;i<l;i++)
- v += alpha[i] * (G[i] + p[i]);
-
- si->obj = v/2;
- }
-
- // put back the solution
- {
- for(int i=0;i<l;i++)
- alpha_[active_set[i]] = alpha[i];
- }
-
- // juggle everything back
- /*{
- for(int i=0;i<l;i++)
- while(active_set[i] != i)
- swap_index(i,active_set[i]);
- // or Q.swap_index(i,active_set[i]);
- }*/
-
- si->upper_bound_p = Cp;
- si->upper_bound_n = Cn;
-
- info("\noptimization finished, #iter = %d\n",iter);
-
- delete[] p;
- delete[] y;
- delete[] alpha;
- delete[] alpha_status;
- delete[] active_set;
- delete[] G;
- delete[] G_bar;
-}
-
-// return 1 if already optimal, return 0 otherwise
-int Solver::select_working_set(int &out_i, int &out_j)
-{
- // return i,j such that
- // i: maximizes -y_i * grad(f)_i, i in I_up(\alpha)
- // j: minimizes the decrease of obj value
- // (if quadratic coefficeint <= 0, replace it with tau)
- // -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha)
-
- double Gmax = -INF;
- double Gmax2 = -INF;
- int Gmax_idx = -1;
- int Gmin_idx = -1;
- double obj_diff_min = INF;
-
- for(int t=0;t<active_size;t++)
- if(y[t]==+1)
- {
- if(!is_upper_bound(t))
- if(-G[t] >= Gmax)
- {
- Gmax = -G[t];
- Gmax_idx = t;
- }
- }
- else
- {
- if(!is_lower_bound(t))
- if(G[t] >= Gmax)
- {
- Gmax = G[t];
- Gmax_idx = t;
- }
- }
-
- int i = Gmax_idx;
- const Qfloat *Q_i = NULL;
- if(i != -1) // NULL Q_i not accessed: Gmax=-INF if i=-1
- Q_i = Q->get_Q(i,active_size);
-
- for(int j=0;j<active_size;j++)
- {
- if(y[j]==+1)
- {
- if (!is_lower_bound(j))
- {
- double grad_diff=Gmax+G[j];
- if (G[j] >= Gmax2)
- Gmax2 = G[j];
- if (grad_diff > 0)
- {
- double obj_diff;
- double quad_coef=Q_i[i]+QD[j]-2*y[i]*Q_i[j];
- if (quad_coef > 0)
- obj_diff = -(grad_diff*grad_diff)/quad_coef;
- else
- obj_diff = -(grad_diff*grad_diff)/TAU;
-
- if (obj_diff <= obj_diff_min)
- {
- Gmin_idx=j;
- obj_diff_min = obj_diff;
- }
- }
- }
- }
- else
- {
- if (!is_upper_bound(j))
- {
- double grad_diff= Gmax-G[j];
- if (-G[j] >= Gmax2)
- Gmax2 = -G[j];
- if (grad_diff > 0)
- {
- double obj_diff;
- double quad_coef=Q_i[i]+QD[j]+2*y[i]*Q_i[j];
- if (quad_coef > 0)
- obj_diff = -(grad_diff*grad_diff)/quad_coef;
- else
- obj_diff = -(grad_diff*grad_diff)/TAU;
-
- if (obj_diff <= obj_diff_min)
- {
- Gmin_idx=j;
- obj_diff_min = obj_diff;
- }
- }
- }
- }
- }
-
- if(Gmax+Gmax2 < eps)
- return 1;
-
- out_i = Gmax_idx;
- out_j = Gmin_idx;
- return 0;
-}
-
-bool Solver::be_shrunken(int i, double Gmax1, double Gmax2)
-{
- if(is_upper_bound(i))
- {
- if(y[i]==+1)
- return(-G[i] > Gmax1);
- else
- return(-G[i] > Gmax2);
- }
- else if(is_lower_bound(i))
- {
- if(y[i]==+1)
- return(G[i] > Gmax2);
- else
- return(G[i] > Gmax1);
- }
- else
- return(false);
-}
-
-void Solver::do_shrinking()
-{
- int i;
- double Gmax1 = -INF; // max { -y_i * grad(f)_i | i in I_up(\alpha) }
- double Gmax2 = -INF; // max { y_i * grad(f)_i | i in I_low(\alpha) }
-
- // find maximal violating pair first
- for(i=0;i<active_size;i++)
- {
- if(y[i]==+1)
- {
- if(!is_upper_bound(i))
- {
- if(-G[i] >= Gmax1)
- Gmax1 = -G[i];
- }
- if(!is_lower_bound(i))
- {
- if(G[i] >= Gmax2)
- Gmax2 = G[i];
- }
- }
- else
- {
- if(!is_upper_bound(i))
- {
- if(-G[i] >= Gmax2)
- Gmax2 = -G[i];
- }
- if(!is_lower_bound(i))
- {
- if(G[i] >= Gmax1)
- Gmax1 = G[i];
- }
- }
- }
-
- // shrink
-
- for(i=0;i<active_size;i++)
- if (be_shrunken(i, Gmax1, Gmax2))
- {
- active_size--;
- while (active_size > i)
- {
- if (!be_shrunken(active_size, Gmax1, Gmax2))
- {
- swap_index(i,active_size);
- break;
- }
- active_size--;
- }
- }
-
- // unshrink, check all variables again before final iterations
-
- if(unshrinked || Gmax1 + Gmax2 > eps*10) return;
-
- unshrinked = true;
- reconstruct_gradient();
-
- for(i=l-1;i>=active_size;i--)
- if (!be_shrunken(i, Gmax1, Gmax2))
- {
- while (active_size < i)
- {
- if (be_shrunken(active_size, Gmax1, Gmax2))
- {
- swap_index(i,active_size);
- break;
- }
- active_size++;
- }
- active_size++;
- }
-}
-
-double Solver::calculate_rho()
-{
- double r;
- int nr_free = 0;
- double ub = INF, lb = -INF, sum_free = 0;
- for(int i=0;i<active_size;i++)
- {
- double yG = y[i]*G[i];
-
- if(is_upper_bound(i))
- {
- if(y[i]==-1)
- ub = min(ub,yG);
- else
- lb = max(lb,yG);
- }
- else if(is_lower_bound(i))
- {
- if(y[i]==+1)
- ub = min(ub,yG);
- else
- lb = max(lb,yG);
- }
- else
- {
- ++nr_free;
- sum_free += yG;
- }
- }
-
- if(nr_free>0)
- r = sum_free/nr_free;
- else
- r = (ub+lb)/2;
-
- return r;
-}
-
-// ****************************************************************************************
-// Solver for nu-svm classification and regression
-// ****************************************************************************************
-// additional constraint: e^T \alpha = constant
-//
-class Solver_NU : public Solver
-{
-public:
- Solver_NU() {}
- void Solve(int l, const QMatrix& Q, const double *p, const schar *y,
- double *alpha, double Cp, double Cn, double eps,
- SolutionInfo* si, int shrinking)
- {
- this->si = si;
- Solver::Solve(l,Q,p,y,alpha,Cp,Cn,eps,si,shrinking);
- }
-private:
- SolutionInfo *si;
- int select_working_set(int &i, int &j);
- double calculate_rho();
- bool be_shrunken(int i, double Gmax1, double Gmax2, double Gmax3, double Gmax4);
- void do_shrinking();
-};
-
-// return 1 if already optimal, return 0 otherwise
-int Solver_NU::select_working_set(int &out_i, int &out_j)
-{
- // return i,j such that y_i = y_j and
- // i: maximizes -y_i * grad(f)_i, i in I_up(\alpha)
- // j: minimizes the decrease of obj value
- // (if quadratic coefficeint <= 0, replace it with tau)
- // -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha)
-
- double Gmaxp = -INF;
- double Gmaxp2 = -INF;
- int Gmaxp_idx = -1;
-
- double Gmaxn = -INF;
- double Gmaxn2 = -INF;
- int Gmaxn_idx = -1;
-
- int Gmin_idx = -1;
- double obj_diff_min = INF;
-
- for(int t=0;t<active_size;t++)
- if(y[t]==+1)
- {
- if(!is_upper_bound(t))
- if(-G[t] >= Gmaxp)
- {
- Gmaxp = -G[t];
- Gmaxp_idx = t;
- }
- }
- else
- {
- if(!is_lower_bound(t))
- if(G[t] >= Gmaxn)
- {
- Gmaxn = G[t];
- Gmaxn_idx = t;
- }
- }
-
- int ip = Gmaxp_idx;
- int in = Gmaxn_idx;
- const Qfloat *Q_ip = NULL;
- const Qfloat *Q_in = NULL;
- if(ip != -1) // NULL Q_ip not accessed: Gmaxp=-INF if ip=-1
- Q_ip = Q->get_Q(ip,active_size);
- if(in != -1)
- Q_in = Q->get_Q(in,active_size);
-
- for(int j=0;j<active_size;j++)
- {
- if(y[j]==+1)
- {
- if (!is_lower_bound(j))
- {
- double grad_diff=Gmaxp+G[j];
- if (G[j] >= Gmaxp2)
- Gmaxp2 = G[j];
- if (grad_diff > 0)
- {
- double obj_diff;
- double quad_coef = Q_ip[ip]+QD[j]-2*Q_ip[j];
- if (quad_coef > 0)
- obj_diff = -(grad_diff*grad_diff)/quad_coef;
- else
- obj_diff = -(grad_diff*grad_diff)/TAU;
-
- if (obj_diff <= obj_diff_min)
- {
- Gmin_idx=j;
- obj_diff_min = obj_diff;
- }
- }
- }
- }
- else
- {
- if (!is_upper_bound(j))
- {
- double grad_diff=Gmaxn-G[j];
- if (-G[j] >= Gmaxn2)
- Gmaxn2 = -G[j];
- if (grad_diff > 0)
- {
- double obj_diff;
- double quad_coef = Q_in[in]+QD[j]-2*Q_in[j];
- if (quad_coef > 0)
- obj_diff = -(grad_diff*grad_diff)/quad_coef;
- else
- obj_diff = -(grad_diff*grad_diff)/TAU;
-
- if (obj_diff <= obj_diff_min)
- {
- Gmin_idx=j;
- obj_diff_min = obj_diff;
- }
- }
- }
- }
- }
-
- if(max(Gmaxp+Gmaxp2,Gmaxn+Gmaxn2) < eps)
- return 1;
-
- if (y[Gmin_idx] == +1)
- out_i = Gmaxp_idx;
- else
- out_i = Gmaxn_idx;
- out_j = Gmin_idx;
-
- return 0;
-}
-
-bool Solver_NU::be_shrunken(int i, double Gmax1, double Gmax2, double Gmax3, double Gmax4)
-{
- if(is_upper_bound(i))
- {
- if(y[i]==+1)
- return(-G[i] > Gmax1);
- else
- return(-G[i] > Gmax4);
- }
- else if(is_lower_bound(i))
- {
- if(y[i]==+1)
- return(G[i] > Gmax2);
- else
- return(G[i] > Gmax3);
- }
- else
- return(false);
-}
-
-void Solver_NU::do_shrinking()
-{
- double Gmax1 = -INF; // max { -y_i * grad(f)_i | y_i = +1, i in I_up(\alpha) }
- double Gmax2 = -INF; // max { y_i * grad(f)_i | y_i = +1, i in I_low(\alpha) }
- double Gmax3 = -INF; // max { -y_i * grad(f)_i | y_i = -1, i in I_up(\alpha) }
- double Gmax4 = -INF; // max { y_i * grad(f)_i | y_i = -1, i in I_low(\alpha) }
-
- // find maximal violating pair first
- int i;
- for(i=0;i<active_size;i++)
- {
- if(!is_upper_bound(i))
- {
- if(y[i]==+1)
- {
- if(-G[i] > Gmax1) Gmax1 = -G[i];
- }
- else if(-G[i] > Gmax4) Gmax4 = -G[i];
- }
- if(!is_lower_bound(i))
- {
- if(y[i]==+1)
- {
- if(G[i] > Gmax2) Gmax2 = G[i];
- }
- else if(G[i] > Gmax3) Gmax3 = G[i];
- }
- }
-
- // shrinking
-
- for(i=0;i<active_size;i++)
- if (be_shrunken(i, Gmax1, Gmax2, Gmax3, Gmax4))
- {
- active_size--;
- while (active_size > i)
- {
- if (!be_shrunken(active_size, Gmax1, Gmax2, Gmax3, Gmax4))
- {
- swap_index(i,active_size);
- break;
- }
- active_size--;
- }
- }
-
- // unshrink, check all variables again before final iterations
-
- if(unshrinked || max(Gmax1+Gmax2,Gmax3+Gmax4) > eps*10) return;
-
- unshrinked = true;
- reconstruct_gradient();
-
- for(i=l-1;i>=active_size;i--)
- if (!be_shrunken(i, Gmax1, Gmax2, Gmax3, Gmax4))
- {
- while (active_size < i)
- {
- if (be_shrunken(active_size, Gmax1, Gmax2, Gmax3, Gmax4))
- {
- swap_index(i,active_size);
- break;
- }
- active_size++;
- }
- active_size++;
- }
-}
-
-double Solver_NU::calculate_rho()
-{
- int nr_free1 = 0,nr_free2 = 0;
- double ub1 = INF, ub2 = INF;
- double lb1 = -INF, lb2 = -INF;
- double sum_free1 = 0, sum_free2 = 0;
-
- for(int i=0;i<active_size;i++)
- {
- if(y[i]==+1)
- {
- if(is_upper_bound(i))
- lb1 = max(lb1,G[i]);
- else if(is_lower_bound(i))
- ub1 = min(ub1,G[i]);
- else
- {
- ++nr_free1;
- sum_free1 += G[i];
- }
- }
- else
- {
- if(is_upper_bound(i))
- lb2 = max(lb2,G[i]);
- else if(is_lower_bound(i))
- ub2 = min(ub2,G[i]);
- else
- {
- ++nr_free2;
- sum_free2 += G[i];
- }
- }
- }
-
- double r1,r2;
- if(nr_free1 > 0)
- r1 = sum_free1/nr_free1;
- else
- r1 = (ub1+lb1)/2;
-
- if(nr_free2 > 0)
- r2 = sum_free2/nr_free2;
- else
- r2 = (ub2+lb2)/2;
-
- si->r = (r1+r2)/2;
- return (r1-r2)/2;
-}
-
-// ****************************************************************************************
-// Q matrices for various formulations
-// ****************************************************************************************
-class SVC_Q: public Kernel
-{
-public:
- SVC_Q(const svm_problem& prob, const svm_parameter& param, const schar *y_)
- :Kernel(prob.l, prob.x, param)
- {
- clone(y,y_,prob.l);
- cache = new Cache(prob.l,(long int)(param.cache_size*(1<<20)));
- QD = new Qfloat[prob.l];
- for(int i=0;i<prob.l;i++)
- QD[i]= (Qfloat)(this->*kernel_function)(i,i,this->m_param);
- }
-
- Qfloat *get_Q(int i, int len) const
- {
- Qfloat *data;
- int start;
- if((start = cache->get_data(i,&data,len)) < len)
- {
- for(int j=start;j<len;j++)
- data[j] = (Qfloat)(y[i]*y[j]*(this->*kernel_function)(i,j,this->m_param));
- }
- return data;
- }
-
- Qfloat *get_QD() const
- {
- return QD;
- }
-
- void swap_index(int i, int j) const
- {
- cache->swap_index(i,j);
- Kernel::swap_index(i,j);
- swap(y[i],y[j]);
- swap(QD[i],QD[j]);
- }
-
- ~SVC_Q()
- {
- delete[] y;
- delete cache;
- delete[] QD;
- }
-private:
- schar *y;
- Cache *cache;
- Qfloat *QD;
-};
-
-
-// ****************************************************************************************
-// ONE_CLASS_Q
-// ****************************************************************************************
-class ONE_CLASS_Q: public Kernel
-{
-public:
- ONE_CLASS_Q(const svm_problem& prob, const svm_parameter& param)
- :Kernel(prob.l, prob.x, param)
- {
- cache = new Cache(prob.l,(long int)(param.cache_size*(1<<20)));
- QD = new Qfloat[prob.l];
- for(int i=0;i<prob.l;i++)
- //OTB's modifications
- QD[i]= (Qfloat)(this->*kernel_function)(i,i,this->m_param);
- }
-
- Qfloat *get_Q(int i, int len) const
- {
- Qfloat *data;
- int start;
- if((start = cache->get_data(i,&data,len)) < len)
- {
- //OTB's modifications
- for(int j=start;j<len;j++)
- data[j] = (Qfloat)(this->*kernel_function)(i,j,this->m_param);
- }
- return data;
- }
-
- Qfloat *get_QD() const
- {
- return QD;
- }
-
- void swap_index(int i, int j) const
- {
- cache->swap_index(i,j);
- Kernel::swap_index(i,j);
- swap(QD[i],QD[j]);
- }
-
- ~ONE_CLASS_Q()
- {
- delete cache;
- delete[] QD;
- }
-private:
- Cache *cache;
- Qfloat *QD;
-};
-
-class SVR_Q: public Kernel
-{
-public:
- SVR_Q(const svm_problem& prob, const svm_parameter& param)
- :Kernel(prob.l, prob.x, param)
- {
- l = prob.l;
- cache = new Cache(l,(long int)(param.cache_size*(1<<20)));
- QD = new Qfloat[2*l];
- sign = new schar[2*l];
- index = new int[2*l];
- for(int k=0;k<l;k++)
- {
- sign[k] = 1;
- sign[k+l] = -1;
- index[k] = k;
- index[k+l] = k;
- //OTB's modifications
- QD[k]= (Qfloat)(this->*kernel_function)(k,k,this->m_param);
- QD[k+l]=QD[k];
- }
- buffer[0] = new Qfloat[2*l];
- buffer[1] = new Qfloat[2*l];
- next_buffer = 0;
- }
-
- void swap_index(int i, int j) const
- {
- swap(sign[i],sign[j]);
- swap(index[i],index[j]);
- swap(QD[i],QD[j]);
- }
-
- Qfloat *get_Q(int i, int len) const
- {
- Qfloat *data;
- int real_i = index[i];
- if(cache->get_data(real_i,&data,l) < l)
- {
- for(int j=0;j<l;j++)
- //OTB's modifications
- data[j] = (Qfloat)(this->*kernel_function)(real_i,j,this->m_param);
- }
-
- // reorder and copy
- Qfloat *buf = buffer[next_buffer];
- next_buffer = 1 - next_buffer;
- schar si = sign[i];
- for(int j=0;j<len;j++)
- buf[j] = si * sign[j] * data[index[j]];
- return buf;
- }
-
- Qfloat *get_QD() const
- {
- return QD;
- }
-
- ~SVR_Q()
- {
- delete cache;
- delete[] sign;
- delete[] index;
- delete[] buffer[0];
- delete[] buffer[1];
- delete[] QD;
- }
-private:
- int l;
- Cache *cache;
- schar *sign;
- int *index;
- mutable int next_buffer;
- Qfloat *buffer[2];
- Qfloat *QD;
-};
-
-// ****************************************************************************************
-// construct and solve various formulations
-// ****************************************************************************************
-static void solve_c_svc(
- const svm_problem *prob, const svm_parameter* param,
- double *alpha, Solver::SolutionInfo* si, double Cp, double Cn)
-{
- int l = prob->l;
- double *minus_ones = new double[l];
- schar *y = new schar[l];
-
- int i;
-
- for(i=0;i<l;i++)
- {
- alpha[i] = 0;
- minus_ones[i] = -1;
- if(prob->y[i] > 0) y[i] = +1; else y[i]=-1;
- }
-
- Solver s;
- s.Solve(l, SVC_Q(*prob,*param,y), minus_ones, y,
- alpha, Cp, Cn, param->eps, si, param->shrinking);
-
- double sum_alpha=0;
- for(i=0;i<l;i++)
- sum_alpha += alpha[i];
-
- if (Cp==Cn)
- info("nu = %f\n", sum_alpha/(Cp*prob->l));
-
- for(i=0;i<l;i++)
- alpha[i] *= y[i];
-
- delete[] minus_ones;
- delete[] y;
-}
-
-static void solve_nu_svc(
- const svm_problem *prob, const svm_parameter *param,
- double *alpha, Solver::SolutionInfo* si)
-{
- int i;
- int l = prob->l;
- double nu = param->nu;
-
- schar *y = new schar[l];
-
- for(i=0;i<l;i++)
- if(prob->y[i]>0)
- y[i] = +1;
- else
- y[i] = -1;
-
- double sum_pos = nu*l/2;
- double sum_neg = nu*l/2;
-
- for(i=0;i<l;i++)
- if(y[i] == +1)
- {
- alpha[i] = min(1.0,sum_pos);
- sum_pos -= alpha[i];
- }
- else
- {
- alpha[i] = min(1.0,sum_neg);
- sum_neg -= alpha[i];
- }
-
- double *zeros = new double[l];
-
- for(i=0;i<l;i++)
- zeros[i] = 0;
-
- Solver_NU s;
- s.Solve(l, SVC_Q(*prob,*param,y), zeros, y,
- alpha, 1.0, 1.0, param->eps, si, param->shrinking);
- double r = si->r;
-
- info("C = %f\n",1/r);
-
- for(i=0;i<l;i++)
- alpha[i] *= y[i]/r;
-
- si->rho /= r;
- si->obj /= (r*r);
- si->upper_bound_p = 1/r;
- si->upper_bound_n = 1/r;
-
- delete[] y;
- delete[] zeros;
-}
-
-// estimatethe support of high-dimensional distribution.
-static void solve_one_class(
- const svm_problem *prob, const svm_parameter *param,
- double *alpha, Solver::SolutionInfo* si)
-{
- int l = prob->l;
- double *zeros = new double[l];
- schar *ones = new schar[l];
- int i;
-
- int n = (int)(param->nu*prob->l); // # of alpha's at upper bound
-
- for(i=0;i<n;i++)
- alpha[i] = 1;
- if(n<prob->l)
- alpha[n] = param->nu * prob->l - n;
- for(i=n+1;i<l;i++)
- alpha[i] = 0;
-
- for(i=0;i<l;i++)
- {
- zeros[i] = 0;
- ones[i] = 1;
- }
-
- Solver s;
- s.Solve(l, ONE_CLASS_Q(*prob,*param), zeros, ones,
- alpha, 1.0, 1.0, param->eps, si, param->shrinking);
-
- delete[] zeros;
- delete[] ones;
-}
-
-
-// epsilon support vector regression
-static void solve_epsilon_svr(
- const svm_problem *prob, const svm_parameter *param,
- double *alpha, Solver::SolutionInfo* si)
-{
- int l = prob->l;
- double *alpha2 = new double[2*l];
- double *linear_term = new double[2*l];
- schar *y = new schar[2*l];
- int i;
-
- for(i=0;i<l;i++)
- {
- alpha2[i] = 0;
- linear_term[i] = param->p - prob->y[i];
- y[i] = 1;
-
- alpha2[i+l] = 0;
- linear_term[i+l] = param->p + prob->y[i];
- y[i+l] = -1;
- }
-
- Solver s;
- s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
- alpha2, param->C, param->C, param->eps, si, param->shrinking);
-
- double sum_alpha = 0;
- for(i=0;i<l;i++)
- {
- alpha[i] = alpha2[i] - alpha2[i+l];
- sum_alpha += fabs(alpha[i]);
- }
- info("nu = %f\n",sum_alpha/(param->C*l));
-
- delete[] alpha2;
- delete[] linear_term;
- delete[] y;
-}
-
-// regression using nu to control the number of support vectors.
-static void solve_nu_svr(
- const svm_problem *prob, const svm_parameter *param,
- double *alpha, Solver::SolutionInfo* si)
-{
- int l = prob->l;
- double C = param->C;
- double *alpha2 = new double[2*l];
- double *linear_term = new double[2*l];
- schar *y = new schar[2*l];
- int i;
-
- double sum = C * param->nu * l / 2;
- for(i=0;i<l;i++)
- {
- alpha2[i] = alpha2[i+l] = min(sum,C);
- sum -= alpha2[i];
-
- linear_term[i] = - prob->y[i];
- y[i] = 1;
-
- linear_term[i+l] = prob->y[i];
- y[i+l] = -1;
- }
-
- Solver_NU s;
- s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
- alpha2, C, C, param->eps, si, param->shrinking);
-
- info("epsilon = %f\n",-si->r);
-
- for(i=0;i<l;i++)
- alpha[i] = alpha2[i] - alpha2[i+l];
-
- delete[] alpha2;
- delete[] linear_term;
- delete[] y;
-}
-
-// ****************************************************************************************
-// decision_function
-// ****************************************************************************************
-struct decision_function
-{
- double *alpha;
- double rho;
-};
-
-decision_function svm_train_one(
- const svm_problem *prob, const svm_parameter *param,
- double Cp, double Cn)
-{
- double *alpha = Malloc(double,prob->l);
- Solver::SolutionInfo si;
- switch(param->svm_type)
- {
- case C_SVC:
- solve_c_svc(prob,param,alpha,&si,Cp,Cn);
- break;
- case NU_SVC:
- solve_nu_svc(prob,param,alpha,&si);
- break;
- case ONE_CLASS:
- solve_one_class(prob,param,alpha,&si);
- break;
- case EPSILON_SVR:
- solve_epsilon_svr(prob,param,alpha,&si);
- break;
- case NU_SVR:
- solve_nu_svr(prob,param,alpha,&si);
- break;
- }
-
- info("obj = %f, rho = %f\n",si.obj,si.rho);
-
- // output SVs
-
- int nSV = 0;
- int nBSV = 0;
- for(int i=0;i<prob->l;i++)
- {
- if(fabs(alpha[i]) > 0)
- {
- ++nSV;
- if(prob->y[i] > 0)
- {
- if(fabs(alpha[i]) >= si.upper_bound_p)
- ++nBSV;
- }
- else
- {
- if(fabs(alpha[i]) >= si.upper_bound_n)
- ++nBSV;
- }
- }
- }
-
- info("nSV = %d, nBSV = %d\n",nSV,nBSV);
-
- decision_function f;
- f.alpha = alpha;
- f.rho = si.rho;
- return f;
-}
-
-
-// Platt's binary SVM Probablistic Output: an improvement from Lin et al.
-void sigmoid_train(
- int l, const double *dec_values, const double *labels,
- double& A, double& B)
-{
- double prior1=0, prior0 = 0;
- int i;
-
- for (i=0;i<l;i++)
- if (labels[i] > 0) prior1+=1;
- else prior0+=1;
-
- int max_iter=100; // Maximal number of iterations
- double min_step=1e-10; // Minimal step taken in line search
- double sigma=1e-3; // For numerically strict PD of Hessian
- double eps=1e-5;
- double hiTarget=(prior1+1.0)/(prior1+2.0);
- double loTarget=1/(prior0+2.0);
- double *t=Malloc(double,l);
- double fApB,p,q,h11,h22,h21,g1,g2,det,dA,dB,gd,stepsize;
- double newA,newB,newf,d1,d2;
- int iter;
-
- // Initial Point and Initial Fun Value
- A=0.0; B=log((prior0+1.0)/(prior1+1.0));
- double fval = 0.0;
-
- for (i=0;i<l;i++)
- {
- if (labels[i]>0) t[i]=hiTarget;
- else t[i]=loTarget;
- fApB = dec_values[i]*A+B;
- if (fApB>=0)
- fval += t[i]*fApB + log(1+exp(-fApB));
- else
- fval += (t[i] - 1)*fApB +log(1+exp(fApB));
- }
- for (iter=0;iter<max_iter;iter++)
- {
- // Update Gradient and Hessian (use H' = H + sigma I)
- h11=sigma; // numerically ensures strict PD
- h22=sigma;
- h21=0.0;g1=0.0;g2=0.0;
- for (i=0;i<l;i++)
- {
- fApB = dec_values[i]*A+B;
- if (fApB >= 0)
- {
- p=exp(-fApB)/(1.0+exp(-fApB));
- q=1.0/(1.0+exp(-fApB));
- }
- else
- {
- p=1.0/(1.0+exp(fApB));
- q=exp(fApB)/(1.0+exp(fApB));
- }
- d2=p*q;
- h11+=dec_values[i]*dec_values[i]*d2;
- h22+=d2;
- h21+=dec_values[i]*d2;
- d1=t[i]-p;
- g1+=dec_values[i]*d1;
- g2+=d1;
- }
-
- // Stopping Criteria
- if (fabs(g1)<eps && fabs(g2)<eps)
- break;
-
- // Finding Newton direction: -inv(H') * g
- det=h11*h22-h21*h21;
- dA=-(h22*g1 - h21 * g2) / det;
- dB=-(-h21*g1+ h11 * g2) / det;
- gd=g1*dA+g2*dB;
-
-
- stepsize = 1; // Line Search
- while (stepsize >= min_step)
- {
- newA = A + stepsize * dA;
- newB = B + stepsize * dB;
-
- // New function value
- newf = 0.0;
- for (i=0;i<l;i++)
- {
- fApB = dec_values[i]*newA+newB;
- if (fApB >= 0)
- newf += t[i]*fApB + log(1+exp(-fApB));
- else
- newf += (t[i] - 1)*fApB +log(1+exp(fApB));
- }
- // Check sufficient decrease
- if (newf<fval+0.0001*stepsize*gd)
- {
- A=newA;B=newB;fval=newf;
- break;
- }
- else
- stepsize = stepsize / 2.0;
- }
-
- if (stepsize < min_step)
- {
- info("Line search fails in two-class probability estimates\n");
- break;
- }
- }
-
- if (iter>=max_iter)
- info("Reaching maximal iterations in two-class probability estimates\n");
- delete [](t);
-}
-
-double sigmoid_predict(double decision_value, double A, double B)
-{
- double fApB = decision_value*A+B;
- if (fApB >= 0)
- return exp(-fApB)/(1.0+exp(-fApB));
- else
- return 1.0/(1+exp(fApB)) ;
-}
-
-// Method 2 from the multiclass_prob paper by Wu, Lin, and Weng
-void multiclass_probability(int k, double **r, double *p)
-{
- int t,j;
- int iter = 0, max_iter=max(100,k);
- double **Q=Malloc(double *,k);
- double *Qp=Malloc(double,k);
- double pQp, eps=0.005/k;
-
- for (t=0;t<k;t++)
- {
- p[t]=1.0/k; // Valid if k = 1
- Q[t]=Malloc(double,k);
- Q[t][t]=0;
- for (j=0;j<t;j++)
- {
- Q[t][t]+=r[j][t]*r[j][t];
- Q[t][j]=Q[j][t];
- }
- for (j=t+1;j<k;j++)
- {
- Q[t][t]+=r[j][t]*r[j][t];
- Q[t][j]=-r[j][t]*r[t][j];
- }
- }
- for (iter=0;iter<max_iter;iter++)
- {
- // stopping condition, recalculate QP,pQP for numerical accuracy
- pQp=0;
- for (t=0;t<k;t++)
- {
- Qp[t]=0;
- for (j=0;j<k;j++)
- Qp[t]+=Q[t][j]*p[j];
- pQp+=p[t]*Qp[t];
- }
- double max_error=0;
- for (t=0;t<k;t++)
- {
- double error=fabs(Qp[t]-pQp);
- if (error>max_error)
- max_error=error;
- }
- if (max_error<eps) break;
-
- for (t=0;t<k;t++)
- {
- double diff=(-Qp[t]+pQp)/Q[t][t];
- p[t]+=diff;
- pQp=(pQp+diff*(diff*Q[t][t]+2*Qp[t]))/(1+diff)/(1+diff);
- for (j=0;j<k;j++)
- {
- Qp[j]=(Qp[j]+diff*Q[t][j])/(1+diff);
- p[j]/=(1+diff);
- }
- }
- }
- if (iter>=max_iter)
- info("Exceeds max_iter in multiclass_prob\n");
- for(t=0;t<k;t++) delete [](Q[t]);
- delete [](Q);
- delete [](Qp);
-}
-
-// Cross-validation decision values for probability estimates
-void svm_binary_svc_probability(
- const svm_problem *prob, const svm_parameter *param,
- double Cp, double Cn, double& probA, double& probB)
-{
- int i;
- int nr_fold = 5;
- int *perm = Malloc(int,prob->l);
- double *dec_values = Malloc(double,prob->l);
-
- // random shuffle
- for(i=0;i<prob->l;i++) perm[i]=i;
- for(i=0;i<prob->l;i++)
- {
- int j = i+rand()%(prob->l-i);
- swap(perm[i],perm[j]);
- }
- for(i=0;i<nr_fold;i++)
- {
- int begin = i*prob->l/nr_fold;
- int end = (i+1)*prob->l/nr_fold;
- int j,k;
- struct svm_problem subprob;
-
- subprob.l = prob->l-(end-begin);
- subprob.x = Malloc(struct svm_node*,subprob.l);
- subprob.y = Malloc(double,subprob.l);
-
- k=0;
- for(j=0;j<begin;j++)
- {
- subprob.x[k] = prob->x[perm[j]];
- subprob.y[k] = prob->y[perm[j]];
- ++k;
- }
- for(j=end;j<prob->l;j++)
- {
- subprob.x[k] = prob->x[perm[j]];
- subprob.y[k] = prob->y[perm[j]];
- ++k;
- }
- int p_count=0,n_count=0;
- for(j=0;j<k;j++)
- if(subprob.y[j]>0)
- p_count++;
- else
- n_count++;
-
- if(p_count==0 && n_count==0)
- for(j=begin;j<end;j++)
- dec_values[perm[j]] = 0;
- else if(p_count > 0 && n_count == 0)
- for(j=begin;j<end;j++)
- dec_values[perm[j]] = 1;
- else if(p_count == 0 && n_count > 0)
- for(j=begin;j<end;j++)
- dec_values[perm[j]] = -1;
- else
- {
- svm_parameter subparam = *param;
- subparam.probability=0;
- subparam.C=1.0;
- subparam.nr_weight=2;
- subparam.weight_label = Malloc(int,2);
- subparam.weight = Malloc(double,2);
- subparam.weight_label[0]=+1;
- subparam.weight_label[1]=-1;
- subparam.weight[0]=Cp;
- subparam.weight[1]=Cn;
- struct svm_model *submodel = svm_train(&subprob,&subparam);
- for(j=begin;j<end;j++)
- {
- svm_predict_values(submodel,prob->x[perm[j]],&(dec_values[perm[j]]));
- // ensure +1 -1 order; reason not using CV subroutine
- dec_values[perm[j]] *= submodel->label[0];
- }
- svm_destroy_model(submodel);
- svm_destroy_param(&subparam);
- }
- delete [](subprob.x);
- delete [](subprob.y);
- }
- sigmoid_train(prob->l,dec_values,prob->y,probA,probB);
- delete [](dec_values);
- delete [](perm);
-}
-
-// Return parameter of a Laplace distribution
-double svm_svr_probability(
- const svm_problem *prob, const svm_parameter *param)
-{
- int i;
- int nr_fold = 5;
- double *ymv = Malloc(double,prob->l);
- double mae = 0;
-
- svm_parameter newparam = *param;
- newparam.probability = 0;
- svm_cross_validation(prob,&newparam,nr_fold,ymv);
- for(i=0;i<prob->l;i++)
- {
- ymv[i]=prob->y[i]-ymv[i];
- mae += fabs(ymv[i]);
- }
- mae /= prob->l;
- double std=sqrt(2*mae*mae);
- int count=0;
- mae=0;
- for(i=0;i<prob->l;i++)
- if (fabs(ymv[i]) > 5*std)
- count=count+1;
- else
- mae+=fabs(ymv[i]);
- mae /= (prob->l-count);
- info("Prob. model for test data: target value = predicted value + z,\nz: Laplace distribution e^(-|z|/sigma)/(2sigma),sigma= %g\n",mae);
- delete [](ymv);
- return mae;
-}
-
-
-// label: label name, start: begin of each class, count: #data of classes, perm: indices to the original data
-// perm, length l, must be allocated before calling this subroutine
-void svm_group_classes(const svm_problem *prob, int *nr_class_ret, int **label_ret, int **start_ret, int **count_ret, int *perm)
-{
- int l = prob->l;
- int max_nr_class = 16;
- int nr_class = 0;
- int *label = Malloc(int,max_nr_class);
- int *count = Malloc(int,max_nr_class);
- int *data_label = Malloc(int,l);
- int i;
-
- for(i=0;i<l;i++)
- {
- int this_label = (int)prob->y[i];
- int j;
- for(j=0;j<nr_class;j++)
- {
- if(this_label == label[j])
- {
- ++count[j];
- break;
- }
- }
- data_label[i] = j;
- if(j == nr_class)
- {
- if(nr_class == max_nr_class)
- {
- max_nr_class *= 2;
- label = (int *)realloc(label,max_nr_class*sizeof(int));
- count = (int *)realloc(count,max_nr_class*sizeof(int));
- }
- label[nr_class] = this_label;
- count[nr_class] = 1;
- ++nr_class;
- }
- }
-
- int *start = Malloc(int,nr_class);
- start[0] = 0;
- for(i=1;i<nr_class;i++)
- start[i] = start[i-1]+count[i-1];
- for(i=0;i<l;i++)
- {
- perm[start[data_label[i]]] = i;
- ++start[data_label[i]];
- }
- start[0] = 0;
- for(i=1;i<nr_class;i++)
- start[i] = start[i-1]+count[i-1];
-
- *nr_class_ret = nr_class;
- *label_ret = label;
- *start_ret = start;
- *count_ret = count;
- delete [](data_label);
-}
-
-// ****************************************************************************************
-// Interface functions
-// ****************************************************************************************
-svm_model *svm_train(const svm_problem *prob, const svm_parameter *param)
-{
- svm_model *model = Malloc(svm_model,1);
- model->param = *param;
- model->free_sv = 0; // XXX
-
- if(param->svm_type == ONE_CLASS ||
- param->svm_type == EPSILON_SVR ||
- param->svm_type == NU_SVR)
- {
- // regression or one-class-svm
- model->nr_class = 2;
- model->label = NULL;
- model->nSV = NULL;
- model->probA = NULL; model->probB = NULL;
- model->sv_coef = Malloc(double *,1);
-
- if(param->probability &&
- (param->svm_type == EPSILON_SVR ||
- param->svm_type == NU_SVR))
- {
- model->probA = Malloc(double,1);
- model->probA[0] = svm_svr_probability(prob,param);
- }
-
- decision_function f = svm_train_one(prob,param,0,0);
- model->rho = Malloc(double,1);
- model->rho[0] = f.rho;
-
- int nSV = 0;
- int i;
- for(i=0;i<prob->l;i++)
- if(fabs(f.alpha[i]) > 0) ++nSV;
- model->l = nSV;
- model->SV = Malloc(svm_node *,nSV);
- model->sv_coef[0] = Malloc(double,nSV);
- int j = 0;
- for(i=0;i<prob->l;i++)
- if(fabs(f.alpha[i]) > 0)
- {
- model->SV[j] = prob->x[i];
- model->sv_coef[0][j] = f.alpha[i];
- ++j;
- }
-
- delete [](f.alpha);
- }
- else
- {
- // classification
- int l = prob->l;
- int nr_class;
- int *label = NULL;
- int *start = NULL;
- int *count = NULL;
- int *perm = Malloc(int,l);
-
- // group training data of the same class
- svm_group_classes(prob,&nr_class,&label,&start,&count,perm);
- svm_node **x = Malloc(svm_node *,l);
- int i;
- for(i=0;i<l;i++)
- x[i] = prob->x[perm[i]];
-
- // calculate weighted C
-
- double *weighted_C = Malloc(double, nr_class);
- for(i=0;i<nr_class;i++)
- weighted_C[i] = param->C;
- for(i=0;i<param->nr_weight;i++)
- {
- int j;
- for(j=0;j<nr_class;j++)
- if(param->weight_label[i] == label[j])
- break;
- if(j == nr_class)
- fprintf(stderr,"warning: class label %d specified in weight is not found\n", param->weight_label[i]);
- else
- weighted_C[j] *= param->weight[i];
- }
-
- // train k*(k-1)/2 models
-
- bool *nonzero = Malloc(bool,l);
- for(i=0;i<l;i++)
- nonzero[i] = false;
- decision_function *f = Malloc(decision_function,nr_class*(nr_class-1)/2);
-
- double *probA=NULL,*probB=NULL;
- if (param->probability)
- {
- probA=Malloc(double,nr_class*(nr_class-1)/2);
- probB=Malloc(double,nr_class*(nr_class-1)/2);
- }
-
- int p = 0;
- for(i=0;i<nr_class;i++)
- for(int j=i+1;j<nr_class;j++)
- {
- svm_problem sub_prob;
- int si = start[i], sj = start[j];
- int ci = count[i], cj = count[j];
- sub_prob.l = ci+cj;
- sub_prob.x = Malloc(svm_node *,sub_prob.l);
- sub_prob.y = Malloc(double,sub_prob.l);
- int k;
- for(k=0;k<ci;k++)
- {
- sub_prob.x[k] = x[si+k];
- sub_prob.y[k] = +1;
- }
- for(k=0;k<cj;k++)
- {
- sub_prob.x[ci+k] = x[sj+k];
- sub_prob.y[ci+k] = -1;
- }
-
- if(param->probability)
- svm_binary_svc_probability(&sub_prob,param,weighted_C[i],weighted_C[j],probA[p],probB[p]);
-
- f[p] = svm_train_one(&sub_prob,param,weighted_C[i],weighted_C[j]);
- for(k=0;k<ci;k++)
- if(!nonzero[si+k] && fabs(f[p].alpha[k]) > 0)
- nonzero[si+k] = true;
- for(k=0;k<cj;k++)
- if(!nonzero[sj+k] && fabs(f[p].alpha[ci+k]) > 0)
- nonzero[sj+k] = true;
- delete [](sub_prob.x);
- delete [](sub_prob.y);
- ++p;
- }
-
- // build output
-
- model->nr_class = nr_class;
-
- model->label = Malloc(int,nr_class);
- for(i=0;i<nr_class;i++)
- model->label[i] = label[i];
-
- model->rho = Malloc(double,nr_class*(nr_class-1)/2);
- for(i=0;i<nr_class*(nr_class-1)/2;i++)
- model->rho[i] = f[i].rho;
-
- if(param->probability)
- {
- model->probA = Malloc(double,nr_class*(nr_class-1)/2);
- model->probB = Malloc(double,nr_class*(nr_class-1)/2);
- for(i=0;i<nr_class*(nr_class-1)/2;i++)
- {
- model->probA[i] = probA[i];
- model->probB[i] = probB[i];
- }
- }
- else
- {
- model->probA=NULL;
- model->probB=NULL;
- }
-
- int total_sv = 0;
- int *nz_count = Malloc(int,nr_class);
- model->nSV = Malloc(int,nr_class);
- for(i=0;i<nr_class;i++)
- {
- int nSV = 0;
- for(int j=0;j<count[i];j++)
- if(nonzero[start[i]+j])
- {
- ++nSV;
- ++total_sv;
- }
- model->nSV[i] = nSV;
- nz_count[i] = nSV;
- }
-
- info("Total nSV = %d\n",total_sv);
-
- model->l = total_sv;
- model->SV = Malloc(svm_node *,total_sv);
- p = 0;
- for(i=0;i<l;i++)
- if(nonzero[i]) model->SV[p++] = x[i];
-
- int *nz_start = Malloc(int,nr_class);
- nz_start[0] = 0;
- for(i=1;i<nr_class;i++)
- nz_start[i] = nz_start[i-1]+nz_count[i-1];
-
- model->sv_coef = Malloc(double *,nr_class-1);
- for(i=0;i<nr_class-1;i++)
- model->sv_coef[i] = Malloc(double,total_sv);
-
- p = 0;
- for(i=0;i<nr_class;i++)
- for(int j=i+1;j<nr_class;j++)
- {
- // classifier (i,j): coefficients with
- // i are in sv_coef[j-1][nz_start[i]...],
- // j are in sv_coef[i][nz_start[j]...]
-
- int si = start[i];
- int sj = start[j];
- int ci = count[i];
- int cj = count[j];
-
- int q = nz_start[i];
- int k;
- for(k=0;k<ci;k++)
- if(nonzero[si+k])
- model->sv_coef[j-1][q++] = f[p].alpha[k];
- q = nz_start[j];
- for(k=0;k<cj;k++)
- if(nonzero[sj+k])
- model->sv_coef[i][q++] = f[p].alpha[ci+k];
- ++p;
- }
-
- delete [](label);
- delete [](probA);
- delete [](probB);
- delete [](count);
- delete [](perm);
- delete [](start);
- delete [](x);
- delete [](weighted_C);
- delete [](nonzero);
- for(i=0;i<nr_class*(nr_class-1)/2;i++)
- delete [](f[i].alpha);
- delete [](f);
- delete [](nz_count);
- delete [](nz_start);
- }
- return model;
-}
-
-// Stratified cross validation
-void svm_cross_validation(const svm_problem *prob, const svm_parameter *param, int nr_fold, double *target)
-{
- int i;
- int *fold_start = Malloc(int,nr_fold+1);
- int l = prob->l;
- int *perm = Malloc(int,l);
- int nr_class;
-
- // stratified cv may not give leasve-one-out rate
- // Each class to l folds -> some folds may have zero elements
- if((param->svm_type == C_SVC ||
- param->svm_type == NU_SVC) && nr_fold < l)
- {
- int *start = NULL;
- int *label = NULL;
- int *count = NULL;
- svm_group_classes(prob,&nr_class,&label,&start,&count,perm);
-
- // random shuffle and then data grouped by fold using the array perm
- int *fold_count = Malloc(int,nr_fold);
- int c;
- int *index = Malloc(int,l);
- for(i=0;i<l;i++)
- index[i]=perm[i];
- for (c=0; c<nr_class; c++)
- for(i=0;i<count[c];i++)
- {
- int j = i+rand()%(count[c]-i);
- swap(index[start[c]+j],index[start[c]+i]);
- }
- for(i=0;i<nr_fold;i++)
- {
- fold_count[i] = 0;
- for (c=0; c<nr_class;c++)
- fold_count[i]+=(i+1)*count[c]/nr_fold-i*count[c]/nr_fold;
- }
- fold_start[0]=0;
- for (i=1;i<=nr_fold;i++)
- fold_start[i] = fold_start[i-1]+fold_count[i-1];
- for (c=0; c<nr_class;c++)
- for(i=0;i<nr_fold;i++)
- {
- int begin = start[c]+i*count[c]/nr_fold;
- int end = start[c]+(i+1)*count[c]/nr_fold;
- for(int j=begin;j<end;j++)
- {
- perm[fold_start[i]] = index[j];
- fold_start[i]++;
- }
- }
- fold_start[0]=0;
- for (i=1;i<=nr_fold;i++)
- fold_start[i] = fold_start[i-1]+fold_count[i-1];
- delete [](start);
- delete [](label);
- delete [](count);
- delete [](index);
- delete [](fold_count);
- }
- else
- {
- for(i=0;i<l;i++) perm[i]=i;
- for(i=0;i<l;i++)
- {
- int j = i+rand()%(l-i);
- swap(perm[i],perm[j]);
- }
- for(i=0;i<=nr_fold;i++)
- fold_start[i]=i*l/nr_fold;
- }
-
- for(i=0;i<nr_fold;i++)
- {
- int begin = fold_start[i];
- int end = fold_start[i+1];
- int j,k;
- struct svm_problem subprob;
-
- subprob.l = l-(end-begin);
- subprob.x = Malloc(struct svm_node*,subprob.l);
- subprob.y = Malloc(double,subprob.l);
-
- k=0;
- for(j=0;j<begin;j++)
- {
- subprob.x[k] = prob->x[perm[j]];
- subprob.y[k] = prob->y[perm[j]];
- ++k;
- }
- for(j=end;j<l;j++)
- {
- subprob.x[k] = prob->x[perm[j]];
- subprob.y[k] = prob->y[perm[j]];
- ++k;
- }
- struct svm_model *submodel = svm_train(&subprob,param);
- if(param->probability &&
- (param->svm_type == C_SVC || param->svm_type == NU_SVC))
- {
- double *prob_estimates=Malloc(double,svm_get_nr_class(submodel));
- for(j=begin;j<end;j++)
- target[perm[j]] = svm_predict_probability(submodel,prob->x[perm[j]],prob_estimates);
- delete [](prob_estimates);
- }
- else
- for(j=begin;j<end;j++)
- target[perm[j]] = svm_predict(submodel,prob->x[perm[j]]);
- svm_destroy_model(submodel);
- delete [](subprob.x);
- delete [](subprob.y);
- }
- delete [](fold_start);
- delete [](perm);
-}
-
-
-int svm_get_svm_type(const svm_model *model)
-{
- return model->param.svm_type;
-}
-
-int svm_get_nr_class(const svm_model *model)
-{
- return model->nr_class;
-}
-
-void svm_get_labels(const svm_model *model, int* label)
-{
- if (model->label != NULL)
- for(int i=0;i<model->nr_class;i++)
- label[i] = model->label[i];
-}
-
-double svm_get_svr_probability(const svm_model *model)
-{
- if ((model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) &&
- model->probA!=NULL)
- return model->probA[0];
- else
- {
- info("Model doesn't contain information for SVR probability inference\n");
- return 0;
- }
-}
-
-void svm_predict_values(const svm_model *model, const svm_node *x, double* dec_values)
-{
- if(model->param.svm_type == ONE_CLASS ||
- model->param.svm_type == EPSILON_SVR ||
- model->param.svm_type == NU_SVR)
- {
- double *sv_coef = model->sv_coef[0];
- double sum = 0;
- for(int i=0;i<model->l;i++)
- sum += sv_coef[i] * Kernel::k_function(x,model->SV[i],model->param);
- sum -= model->rho[0];
- *dec_values = sum;
- }
- else
- {
- int i;
- int nr_class = model->nr_class;
- int l = model->l;
-
- double *kvalue = Malloc(double,l);
- for(i=0;i<l;i++)
- kvalue[i] = Kernel::k_function(x,model->SV[i],model->param);
-
- int *start = Malloc(int,nr_class);
- start[0] = 0;
- for(i=1;i<nr_class;i++)
- start[i] = start[i-1]+model->nSV[i-1];
-
- int p=0;
- for(i=0;i<nr_class;i++)
- for(int j=i+1;j<nr_class;j++)
- {
- double sum = 0;
- int si = start[i];
- int sj = start[j];
- int ci = model->nSV[i];
- int cj = model->nSV[j];
-
- int k;
- double *coef1 = model->sv_coef[j-1];
- double *coef2 = model->sv_coef[i];
- for(k=0;k<ci;k++)
- sum += coef1[si+k] * kvalue[si+k];
- for(k=0;k<cj;k++)
- sum += coef2[sj+k] * kvalue[sj+k];
- sum -= model->rho[p];
- dec_values[p] = sum;
- p++;
- }
-
- delete [](kvalue);
- delete [](start);
- }
-}
-
-double svm_predict(const svm_model *model, const svm_node *x)
-{
- if(model->param.svm_type == ONE_CLASS ||
- model->param.svm_type == EPSILON_SVR ||
- model->param.svm_type == NU_SVR)
- {
- double res;
- svm_predict_values(model, x, &res);
-
- if(model->param.svm_type == ONE_CLASS)
- return (res>0)?1:-1;
- else
- return res;
- }
- else
- {
- int i;
- int nr_class = model->nr_class;
- double *dec_values = Malloc(double, nr_class*(nr_class-1)/2);
- svm_predict_values(model, x, dec_values);
-
- int *vote = Malloc(int,nr_class);
- for(i=0;i<nr_class;i++)
- vote[i] = 0;
- int pos=0;
- for(i=0;i<nr_class;i++)
- for(int j=i+1;j<nr_class;j++)
- {
- if(dec_values[pos++] > 0)
- ++vote[i];
- else
- ++vote[j];
- }
-
- int vote_max_idx = 0;
- for(i=1;i<nr_class;i++)
- if(vote[i] > vote[vote_max_idx])
- vote_max_idx = i;
- delete [](vote);
- delete [](dec_values);
-
- return model->label[vote_max_idx];
- }
-}
-
-double svm_predict_probability(
- const svm_model *model, const svm_node *x, double *prob_estimates)
-{
- if ((model->param.svm_type == C_SVC || model->param.svm_type == NU_SVC) &&
- model->probA!=NULL && model->probB!=NULL)
- {
- int i;
- int nr_class = model->nr_class;
- double *dec_values = Malloc(double, nr_class*(nr_class-1)/2);
- svm_predict_values(model, x, dec_values);
-
- double min_prob=1e-7;
- double **pairwise_prob=Malloc(double *,nr_class);
- for(i=0;i<nr_class;i++)
- pairwise_prob[i]=Malloc(double,nr_class);
- int k=0;
- for(i=0;i<nr_class;i++)
- for(int j=i+1;j<nr_class;j++)
- {
- pairwise_prob[i][j]=min(max(sigmoid_predict(dec_values[k],model->probA[k],model->probB[k]),min_prob),1-min_prob);
- pairwise_prob[j][i]=1-pairwise_prob[i][j];
- k++;
- }
- multiclass_probability(nr_class,pairwise_prob,prob_estimates);
-
- int prob_max_idx = 0;
- for(i=1;i<nr_class;i++)
- if(prob_estimates[i] > prob_estimates[prob_max_idx])
- prob_max_idx = i;
- for(i=0;i<nr_class;i++)
- delete [](pairwise_prob[i]);
- delete [](dec_values);
- delete [](pairwise_prob);
- return model->label[prob_max_idx];
- }
- else
- return svm_predict(model, x);
-}
-
-const char *svm_type_table[] =
- {
- "c_svc","nu_svc","one_class","epsilon_svr","nu_svr",NULL
- };
-
-const char *kernel_type_table[]=
- {
- //OTB's modifications
- "linear","polynomial","rbf","sigmoid","precomputed","generic","composed",NULL
- };
-
-int svm_save_model(const char *model_file_name, const svm_model *model)
-{
- FILE *fp = fopen(model_file_name,"w");
- if(fp==NULL) return -1;
-
- const svm_parameter& param = model->param;
-
- fprintf(fp,"svm_type %s\n", svm_type_table[param.svm_type]);
- fprintf(fp,"kernel_type %s\n", kernel_type_table[param.kernel_type]);
-
- //OTB's modifications
- if( param.kernel_type == GENERIC )
- {
- if( param.kernel_generic == NULL )
- {
- fprintf(stderr,"generic kernel functor is not initialized\n");
- return -1;
- }
- //Load generic parameters
- int cr = param.kernel_generic->save_parameters(&fp,"generic_kernel_parameters");
- if( cr != 0 )
- {
- fprintf(stderr,"error while saving generic kernel parameters to the file %s.\n",model_file_name);
- }
- }
- if( param.kernel_type == COMPOSED )
- {
- if( param.kernel_composed == NULL )
- {
- fprintf(stderr,"composed kernel functor is not initialized\n");
- return -1;
- }
- //Load generic parameters
- int cr = param.kernel_composed->save_parameters(&fp,"composed_kernel_parameters");
- if( cr != 0 )
- {
- fprintf(stderr,"error while saving composed kernel parameters to the file %s.\n",model_file_name);
- }
- }
-
- if(param.kernel_type == POLY)
- fprintf(fp,"degree %d\n", param.degree);
-
- if(param.kernel_type == POLY || param.kernel_type == RBF || param.kernel_type == SIGMOID)
- fprintf(fp,"gamma %g\n", param.gamma);
-
- if(param.kernel_type == POLY || param.kernel_type == SIGMOID)
- fprintf(fp,"coef0 %g\n", param.coef0);
-
- int nr_class = model->nr_class;
- int l = model->l;
- fprintf(fp, "nr_class %d\n", nr_class);
- fprintf(fp, "total_sv %d\n",l);
-
- if(model->rho)
- {
- fprintf(fp, "rho");
- for(int i=0;i<nr_class*(nr_class-1)/2;i++)
- fprintf(fp," %.10g ",model->rho[i]);
- fprintf(fp, "\n");
- }
-
- if(model->label)
- {
- fprintf(fp, "label");
- for(int i=0;i<nr_class;i++)
- fprintf(fp," %d",model->label[i]);
- fprintf(fp, "\n");
- }
-
- if(model->probA) // regression has probA only
- {
- fprintf(fp, "probA");
- for(int i=0;i<nr_class*(nr_class-1)/2;i++)
- fprintf(fp," %.10g ",model->probA[i]);
- fprintf(fp, "\n");
- }
- if(model->probB)
- {
- fprintf(fp, "probB");
- for(int i=0;i<nr_class*(nr_class-1)/2;i++)
- fprintf(fp," %.10g ",model->probB[i]);
- fprintf(fp, "\n");
- }
-
- if(model->nSV)
- {
- fprintf(fp, "nr_sv");
- for(int i=0;i<nr_class;i++)
- fprintf(fp," %d",model->nSV[i]);
- fprintf(fp, "\n");
- }
-
- fprintf(fp, "SV\n");
- const double * const *sv_coef = model->sv_coef;
- const svm_node * const *SV = model->SV;
-
- for(int i=0;i<l;i++)
- {
- for(int j=0;j<nr_class-1;j++)
- fprintf(fp, "%.16g ",sv_coef[j][i]);
-
- const svm_node *p = SV[i];
-
- if(param.kernel_type == PRECOMPUTED)
- fprintf(fp,"0:%d ",(int)(p->value));
- else
- while(p->index != -1)
- {
- fprintf(fp,"%d:%.8g ",p->index,p->value);
- p++;
- }
- fprintf(fp, "\n");
- }
- if (ferror(fp) != 0 || fclose(fp) != 0) return -1;
- else return 0;
-}
-
-svm_model *svm_load_model(const char *model_file_name, /*otb::*/GenericKernelFunctorBase * generic_kernel_functor)
-{
- FILE *fp = fopen(model_file_name,"rb");
- if(fp==NULL) return NULL;
-
- // read parameters
-
- svm_model *model = Malloc(svm_model,1);
- svm_parameter& param = model->param;
- model->rho = NULL;
- model->probA = NULL;
- model->probB = NULL;
- model->label = NULL;
- model->nSV = NULL;
- model->delete_composed = false;
-
- char cmd[81];
- while(1)
- {
- fscanf(fp,"%80s",cmd);
-
- if(strcmp(cmd,"svm_type")==0)
- {
- fscanf(fp,"%80s",cmd);
- int i;
- for(i=0;svm_type_table[i];i++)
- {
- if(strcmp(svm_type_table[i],cmd)==0)
- {
- param.svm_type=i;
- break;
- }
- }
- if(svm_type_table[i] == NULL)
- {
- fprintf(stderr,"unknown svm type.\n");
- delete [](model->rho);
- delete [](model->label);
- delete [](model->nSV);
- delete [](model);
- return NULL;
- }
- }
- else if(strcmp(cmd,"kernel_type")==0)
- {
- fscanf(fp,"%80s",cmd);
- int i;
- for(i=0;kernel_type_table[i];i++)
- {
- if(strcmp(kernel_type_table[i],cmd)==0)
- {
- param.kernel_type=i;
- break;
- }
- }
- if(kernel_type_table[i] == NULL)
- {
- fprintf(stderr,"unknown kernel function.\n");
- delete [](model->rho);
- delete [](model->label);
- delete [](model->nSV);
- delete [](model);
- return NULL;
- }
- }
- else if(strcmp(cmd,"degree")==0)
- fscanf(fp,"%d",¶m.degree);
- else if(strcmp(cmd,"gamma")==0)
- fscanf(fp,"%lf",¶m.gamma);
- else if(strcmp(cmd,"coef0")==0)
- fscanf(fp,"%lf",¶m.coef0);
- else if(strcmp(cmd,"nr_class")==0)
- fscanf(fp,"%d",&model->nr_class);
- else if(strcmp(cmd,"total_sv")==0)
- fscanf(fp,"%d",&model->l);
- else if(strcmp(cmd,"rho")==0)
- {
- int n = model->nr_class * (model->nr_class-1)/2;
- model->rho = Malloc(double,n);
- for(int i=0;i<n;i++)
- fscanf(fp,"%lf",&model->rho[i]);
- }
- else if(strcmp(cmd,"label")==0)
- {
- int n = model->nr_class;
- model->label = Malloc(int,n);
- for(int i=0;i<n;i++)
- fscanf(fp,"%d",&model->label[i]);
- }
- else if(strcmp(cmd,"probA")==0)
- {
- int n = model->nr_class * (model->nr_class-1)/2;
- model->probA = Malloc(double,n);
- for(int i=0;i<n;i++)
- fscanf(fp,"%lf",&model->probA[i]);
- }
- else if(strcmp(cmd,"probB")==0)
- {
- int n = model->nr_class * (model->nr_class-1)/2;
- model->probB = Malloc(double,n);
- for(int i=0;i<n;i++)
- fscanf(fp,"%lf",&model->probB[i]);
- }
- else if(strcmp(cmd,"nr_sv")==0)
- {
- int n = model->nr_class;
- model->nSV = Malloc(int,n);
- for(int i=0;i<n;i++)
- fscanf(fp,"%d",&model->nSV[i]);
- }
- else if(strcmp(cmd,"SV")==0)
- {
- while(1)
- {
- int c = getc(fp);
- if(c==EOF || c=='\n') break;
- }
- break;
- }
- //OTB's modifications
-
- else if(strcmp(cmd,"generic_kernel_parameters")==0)
- {
- if( param.kernel_type == GENERIC )
- {
- if( generic_kernel_functor == NULL )
- {
- fprintf(stderr,"generic kernel functor is not initialized\n",cmd);
- return NULL;
- }
- param.kernel_generic = generic_kernel_functor;
- //Load generic parameters
- int cr = param.kernel_generic->load_parameters(&fp);
- if( cr != 0 )
- {
- fprintf(stderr,"error while loading generic kernel parameters from the file %s.\n",model_file_name);
- }
- }
- else
- {
- //Read the generic_kernel_parameters line
- fgets(cmd,80,fp);
- }
- }
- else if(strcmp(cmd,"composed_kernel_parameters")==0)
- {
- if( param.kernel_type == COMPOSED )
- {
- //Load generic parameters
- delete generic_kernel_functor;
-
- param.kernel_composed = new ComposedKernelFunctor;
- int cr = param.kernel_composed->load_parameters(&fp);
- model->delete_composed = true;
- //int cr = param.kernel_generic->load_parameters(&fp);
- if( cr != 0 )
- {
- fprintf(stderr,"error while loading composed kernel parameters from the file %s.\n",model_file_name);
- }
- }
- else
- {
- //Read the composed_kernel_parameters line
- fgets(cmd,80,fp);
- }
- }
- else
- {
- fprintf(stderr,"unknown text in model file: [%s]\n",cmd);
- delete [](model->rho);
- delete [](model->label);
- delete [](model->nSV);
- delete [](model);
- return NULL;
- }
- }
- // read sv_coef and SV
-
- int elements = 0;
- long pos = ftell(fp);
-
- while(1)
- {
- int c = fgetc(fp);
- switch(c)
- {
- case '\n':
- // count the '-1' element
- case ':':
- ++elements;
- break;
- case EOF:
- goto out;
- default:
- ;
- }
- }
- out:
- fseek(fp,pos,SEEK_SET);
-
- int m = model->nr_class - 1;
- int l = model->l;
- model->sv_coef = Malloc(double *,m);
- int i;
- for(i=0;i<m;i++)
- model->sv_coef[i] = Malloc(double,l);
- model->SV = Malloc(svm_node*,l);
- for(int n = 0;n<l;++n)
- {
- model->SV[n]=Malloc(svm_node,1);
- model->SV[n]->index = -1;
- model->SV[n]->value = 0.;
- }
- svm_node *x_space=NULL;
- if(l>0)
- {
- x_space = Malloc(svm_node,elements);
- for(int n = 0;n<elements;++n)
- {
- x_space[n].index = -1;
- x_space[n].value = 0.;
- }
- }
-
- int j=0;
- for(i=0;i<l;i++)
- {
- model->SV[i] = &x_space[j];
- for(int k=0;k<m;k++)
- fscanf(fp,"%lf",&model->sv_coef[k][i]);
- while(1)
- {
- int c;
- do {
- c = getc(fp);
- if(c=='\n') goto out2;
- } while(isspace(c));
- ungetc(c,fp);
- fscanf(fp,"%d:%lf",&(x_space[j].index),&(x_space[j].value));
- ++j;
- }
- out2:
- x_space[j++].index = -1;
- }
- if (ferror(fp) != 0 || fclose(fp) != 0) return NULL;
-
- model->free_sv = 1; // XXX
- return model;
-}
-
-//************************************************//
-// OTB's modifications : fonction entiere ajoutee //
-//************************************************//
-svm_model *svm_copy_model( const svm_model *model )
-{
- const svm_parameter& param = model->param;
-
- // instanciated the copy
- svm_model *modelCpy = Malloc(svm_model,1);
- svm_parameter& paramCpy = modelCpy->param;
- modelCpy->rho = NULL;
- modelCpy->probA = NULL;
- modelCpy->probB = NULL;
- modelCpy->label = NULL;
- modelCpy->nSV = NULL;
- modelCpy->delete_composed = false;
-
- // SVM type copy
- paramCpy.svm_type = param.svm_type;
- // Kernel type copy
- paramCpy.kernel_type = param.kernel_type;
- // Param copy
- paramCpy.degree = param.degree;
- paramCpy.gamma = param.gamma;
- paramCpy.coef0 = param.coef0;
- // Model variable
- int nr_class = model->nr_class;
- int l = model->l;
-
- modelCpy->nr_class = nr_class;
- modelCpy->l = l;
- if(model->rho)
- {
- int n = model->nr_class * (model->nr_class-1)/2;
- modelCpy->rho = Malloc(double,n);
- for(int i=0; i<n; i++)
- modelCpy->rho[i] = model->rho[i];
- }
- if(model->label)
- {
- modelCpy->label = Malloc(int,nr_class);
- for(int i=0; i<nr_class; i++)
- modelCpy->label[i] = model->label[i];
- }
- if(model->probA)
- {
- int n = nr_class * (nr_class-1)/2;
- modelCpy->probA = Malloc(double,n);
- for(int i=0; i<n; i++)
- modelCpy->probA[i] = model->probA[i];
- }
- if(model->probB)
- {
- int n = nr_class * (nr_class-1)/2;
- modelCpy->probB = Malloc(double,n);
- for(int i=0; i<n; i++)
- modelCpy->probB[i] = model->probB[i];
- }
- if(model->nSV)
- {
- modelCpy->nSV = Malloc(int,nr_class);
- for(int i=0;i<nr_class;i++)
- modelCpy->nSV[i] = model->nSV[i];
- }
-
- // SV copy
- const double * const *sv_coef = model->sv_coef;
- const svm_node * const *SV = model->SV;
-
-
- modelCpy->SV = Malloc(svm_node*,l);
- svm_node **SVCpy = modelCpy->SV;
-
- modelCpy->sv_coef = Malloc(double *,nr_class-1);
-
- for(int i=0; i<nr_class-1; i++)
- modelCpy->sv_coef[i] = Malloc(double,l);
-
-
- // Compute the total number of SV elements.
- unsigned int elements = 0;
- for (int p=0; p<l; p++)
- {
- const svm_node *tempNode = SV[p];
- while(tempNode->index != -1)
- {
- tempNode++;
- elements++;
- }
- elements++;// for -1 values
- }
-
- if(l>0)
- {
- modelCpy->SV[0] = Malloc(svm_node,elements);
- memcpy( modelCpy->SV[0],model->SV[0],sizeof(svm_node*)*elements);
- }
- svm_node *x_space = modelCpy->SV[0];
-
-
- int j = 0;
- for(int i=0; i<l; i++)
- {
- // sv_coef
- for(int k=0; k<nr_class-1; k++)
- modelCpy->sv_coef[k][i] = sv_coef[k][i];
-
- // SV
- modelCpy->SV[i] = &x_space[j];
- const svm_node *p = SV[i];
- svm_node *pCpy = SVCpy[i];
- while(p->index != -1)
- {
- pCpy->index = p->index;
- pCpy->value = p->value;
- p++;
- pCpy++;
- j++;
- }
- pCpy->index = -1;
- j++;
- }
-
- // Generic kernel copy
- if( param.kernel_type == GENERIC )
- {
- paramCpy.kernel_generic = Malloc(GenericKernelFunctorBase, sizeof(*(param.kernel_generic)));
- memcpy(paramCpy.kernel_generic,param.kernel_generic,sizeof(*(param.kernel_generic)));
- }
- // Composrd kernel copy
- if( param.kernel_type == COMPOSED )
- {
- paramCpy.kernel_composed = Malloc(ComposedKernelFunctor, 1);
- *(paramCpy.kernel_composed) = *(param.kernel_composed);
- }
-
- return modelCpy;
-}
-
-
-
-void svm_destroy_model(svm_model* model)
-{
- if(model->free_sv && model->l > 0)
- delete []((void *)(model->SV[0]));
- for(int i=0;i<model->nr_class-1;i++)
- delete [](model->sv_coef[i]);
- if ( model->delete_composed == true)
- delete model->param.kernel_composed;
- delete [](model->SV);
- delete [](model->sv_coef);
- delete [](model->rho);
- delete [](model->label);
- delete [](model->probA);
- delete [](model->probB);
- delete [](model->nSV);
- delete [](model);
-}
-
-void svm_destroy_param(svm_parameter* param)
-{
- delete [](param->weight_label);
- delete [](param->weight);
-}
-
-const char *svm_check_parameter(const svm_problem *prob, const svm_parameter *param)
-{
- // svm_type
-
- int svm_type = param->svm_type;
- if(svm_type != C_SVC &&
- svm_type != NU_SVC &&
- svm_type != ONE_CLASS &&
- svm_type != EPSILON_SVR &&
- svm_type != NU_SVR)
- return "unknown svm type";
-
- // kernel_type, degree
-
- int kernel_type = param->kernel_type;
- if(kernel_type != LINEAR &&
- kernel_type != POLY &&
- kernel_type != RBF &&
- kernel_type != SIGMOID &&
- kernel_type != PRECOMPUTED &&
- //OTB's modifications
- kernel_type != GENERIC &&
- kernel_type != COMPOSED)
- return "unknown kernel type";
- if ( kernel_type == GENERIC )
- {
- if( param->kernel_generic == NULL )
- return "Generic kernel functor not initialized";
- }
- if ( kernel_type == COMPOSED )
- {
- if( param->kernel_composed == NULL )
- return "Composed kernel functor not initialized";
- }
-
-
-
- if(param->degree < 0)
- return "degree of polynomial kernel < 0";
-
- // cache_size,eps,C,nu,p,shrinking
-
- if(param->cache_size <= 0)
- return "cache_size <= 0";
-
- if(param->eps <= 0)
- return "eps <= 0";
-
- if(svm_type == C_SVC ||
- svm_type == EPSILON_SVR ||
- svm_type == NU_SVR)
- if(param->C <= 0)
- return "C <= 0";
-
- if(svm_type == NU_SVC ||
- svm_type == ONE_CLASS ||
- svm_type == NU_SVR)
- if(param->nu <= 0 || param->nu > 1)
- return "nu <= 0 or nu > 1";
-
- if(svm_type == EPSILON_SVR)
- if(param->p < 0)
- return "p < 0";
-
- if(param->shrinking != 0 &&
- param->shrinking != 1)
- return "shrinking != 0 and shrinking != 1";
-
- if(param->probability != 0 &&
- param->probability != 1)
- return "probability != 0 and probability != 1";
-
- if(param->probability == 1 &&
- svm_type == ONE_CLASS)
- return "one-class SVM probability output not supported yet";
-
-
- // check whether nu-svc is feasible
-
- if(svm_type == NU_SVC)
- {
- int l = prob->l;
- int max_nr_class = 16;
- int nr_class = 0;
- int *label = Malloc(int,max_nr_class);
- int *count = Malloc(int,max_nr_class);
-
- int i;
- for(i=0;i<l;i++)
- {
- int this_label = (int)prob->y[i];
- int j;
- for(j=0;j<nr_class;j++)
- if(this_label == label[j])
- {
- ++count[j];
- break;
- }
- if(j == nr_class)
- {
- if(nr_class == max_nr_class)
- {
- max_nr_class *= 2;
- label = (int *)realloc(label,max_nr_class*sizeof(int));
- count = (int *)realloc(count,max_nr_class*sizeof(int));
- }
- label[nr_class] = this_label;
- count[nr_class] = 1;
- ++nr_class;
- }
- }
-
- for(i=0;i<nr_class;i++)
- {
- int n1 = count[i];
- for(int j=i+1;j<nr_class;j++)
- {
- int n2 = count[j];
- if(param->nu*(n1+n2)/2 > min(n1,n2))
- {
- delete [](label);
- delete [](count);
- return "specified nu is infeasible";
- }
- }
- }
- delete [](label);
- delete [](count);
- }
-
- return NULL;
-}
-
-int svm_check_probability_model(const svm_model *model)
-{
- return ((model->param.svm_type == C_SVC || model->param.svm_type == NU_SVC) &&
- model->probA!=NULL && model->probB!=NULL) ||
- ((model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) &&
- model->probA!=NULL);
-}
-
-
-//OTB's modifications
-// Methods of GenericKernelFunctorBase class
-/*
- namespace otb
- {
-*/
-
-GenericKernelFunctorBase::GenericKernelFunctorBase(const GenericKernelFunctorBase& copy)
-{
- *this = copy;
-}
-
-GenericKernelFunctorBase&
-GenericKernelFunctorBase::operator=(const GenericKernelFunctorBase& copy)
-{
- this->m_MapParameters = copy.m_MapParameters;
- this->m_Name = copy.m_Name;
- return *this;
-}
-
-int
-GenericKernelFunctorBase::
-load_parameters(FILE ** pfile)
-{
- int NbParams(0);
- char keyword[81];
- char value[81];
- // Read functor name
- fscanf(*pfile,"%80s",keyword);
- m_Name = std::string(keyword);
- // Read number of parameters
- fscanf(*pfile,"%d",&NbParams);
- // if( NbParams == 0 ) return -1;
- for ( int cpt=0 ; cpt < NbParams ; cpt++)
- {
- fscanf(*pfile,"%80s",keyword);
- fscanf(*pfile,"%80s",value);
- m_MapParameters[std::string(keyword)] = std::string(value);
- }
- return 0;
-}
-int
-GenericKernelFunctorBase::
-save_parameters(FILE ** pfile, const char * generic_kernel_parameters_keyword)const
-{
- MapConstIterator iter=m_MapParameters.begin();
-
- std::string line(generic_kernel_parameters_keyword);
- std::string strNbParams;
- ::otb::StringStream flux;
- flux << m_MapParameters.size();
- flux >> strNbParams;
- line = line + " " + m_Name + " " + strNbParams;
-
- // DON'T USE ITER because for a COPY ONLY, whereas map.size() IS GOOD, map.end()
- // doesn't give the good answer -> SEGFAULT
- // while( iter != m_MapParameters.end() )
- // {
- // line = line + " " + iter->first + " " + iter->second;
- // ++iter;
- // }
- for ( unsigned int i = 0; i<m_MapParameters.size(); i++)
- {
- line = line + " " + iter->first + " " + iter->second;
- ++iter;
- }
- fprintf(*pfile,"%s\n", line.c_str());
- return 0;
-}
-
-void
-GenericKernelFunctorBase::
-print_parameters(void)const
-{
- MapConstIterator iter=this->GetMapParameters().begin();
-
- std::cout << "Print generic kernel parameters: "<<this->GetMapParameters().size()<<std::endl;
- while( iter != this->GetMapParameters().end() )
- {
- std::cout << " "<<iter->first <<" "<<iter->second<<std::endl;
- ++iter;
- }
-}
-
-double
-GenericKernelFunctorBase::
-dot(const svm_node *px, const svm_node *py)const
-{
- double sum = 0.;
-
- while(px->index != -1 && py->index != -1)
- {
- if(px->index == py->index)
- {
- sum += px->value * py->value;
- ++px;
- ++py;
- }
- else
- {
- if(px->index > py->index)
- ++py;
- else
- ++px;
- }
- }
- return sum;
-}
-
-svm_node *
-GenericKernelFunctorBase::
-sub(const svm_node *px, const svm_node *py) const
- /* compute the difference a-b of two sparse vectors */
- /* Note: SVECTOR lists are not followed, but only the first
- SVECTOR is used */
-{
- long veclength = 1;
-
- const svm_node * pxbis = px;
- const svm_node * pybis = py;
-
- while (px->index != -1 && py->index != -1)
- {
- ++veclength;
- if(px->index == py->index)
- {
- ++px;
- ++py;
- }
- else
- {
- if (px->index < py->index)
- {
- ++px;
- }
- else
- {
- ++py;
- }
- }
- }
-
- while (py->index != -1)
- {
- ++veclength;
- ++py;
- }
- while (px->index != -1)
- {
- ++veclength;
- ++px;
- }
-
- svm_node *vec;
- vec = new svm_node[veclength];
- unsigned long int vecIt = 0;
-
- px = pxbis;
- py = pybis;
-
- while (px->index != -1 && py->index != -1)
- {
- if(px->index == py->index)
- {
- (vec[vecIt])=(*px);
- vec[vecIt].value-=py->value;
-
- if(vec[vecIt].value != 0)
- {
- ++vecIt;
- }
- ++px;
- ++py;
-
- }
- else
- {
- if(px->index > py->index)
- {
- (vec[vecIt])=(*py);
- vec[vecIt].value*=(-1);
- ++vecIt;
- ++py;
- }
- else
- {
- (vec[vecIt])=(*px);
- ++vecIt;
- ++px;
- }
- }
- }
-
- while (py->index != -1)
- {
- (vec[vecIt])=(*py);
- vec[vecIt].value*=(-1);
- ++vecIt;
- ++py;
- }
- while (px->index != -1)
- {
- (vec[vecIt])=(*px);
- ++vecIt;
- ++px;
- }
- return(vec);
-}
-
-
-svm_node *
-GenericKernelFunctorBase::
-add(const svm_node *px, const svm_node *py) const
- /* compute the sum a+b of two sparse vectors */
- /* Note: SVECTOR lists are not followed, but only the first
- SVECTOR is used */
-{
- const svm_node * pxbis = px;
- const svm_node * pybis = py;
- long veclength = 1;
-
- while (px->index != -1 && py->index != -1)
- {
- ++veclength;
- if(px->index == py->index)
- {
- ++px;
- ++py;
- }
- else
- {
- if (px->index < py->index)
- {
- ++px;
- }
- else
- {
- ++py;
- }
- }
- }
-
- while (py->index != -1)
- {
- ++veclength;
- ++py;
- }
- while (px->index != -1)
- {
- ++veclength;
- ++px;
- }
-
- svm_node *vec;
- vec = new svm_node[veclength];
- unsigned long int vecIt = 0;
-
- px = pxbis;
- py = pybis;
-
- while (px->index != -1 && py->index != -1)
- {
- if(px->index == py->index)
- {
- (vec[vecIt])=(*px);
- vec[vecIt].value+=py->value;
- if(vec[vecIt].value != 0)
- {
- ++vecIt;
- }
- ++px;
- ++py;
- }
- else
- {
- if (px->index < py->index)
- {
- (vec[vecIt])=(*px);
- ++vecIt;
- ++px;
- }
- else
- {
- (vec[vecIt])=(*py);
- ++vecIt;
- ++py;
-
- }
- }
- }
- while (py->index != -1)
- {
- (vec[vecIt])=(*py);
- ++vecIt;
- ++py;
- }
- while (px->index != -1)
- {
- (vec[vecIt])=(*px);
- ++vecIt;
- ++px;
- }
-
- return(vec);
-}
-
-// ****************************************************************************************
-// ************************ ComposedKernelFunctor methods ********************/
-// ****************************************************************************************
-
-ComposedKernelFunctor::ComposedKernelFunctor(const ComposedKernelFunctor& copy)
-{
- *this = copy;
-}
-
-ComposedKernelFunctor&
-ComposedKernelFunctor::operator=(const ComposedKernelFunctor& copy)
-{
- // Call Superclass::operator=
- static_cast<Superclass&>(*this) = static_cast<const Superclass&>(copy);
-
- // Copy Self attributes
- this->m_KernelFunctorList = copy.m_KernelFunctorList;
- this->m_HaveToBeDeletedList = copy.m_HaveToBeDeletedList;
- this->m_PonderationList = copy.m_PonderationList;
- return *this;
-}
-
-void
-ComposedKernelFunctor
-::print_parameters(void)const
-{
- MapConstIterator iter = this->GetMapParameters().begin();
- std::cout << "Print composed kernel parameters: "<<this->GetName()<<", "<<this->GetMapParameters().size()<<std::endl;
- while( iter != this->GetMapParameters().end() )
- {
- std::cout << " "<<iter->first <<" "<<iter->second<<std::endl;
- ++iter;
- }
- std::cout<<std::endl;
- std::cout<<"Composition kernels:"<<std::endl;
- if (m_KernelFunctorList.size() != 0 && m_PonderationList.size() != 0 && m_KernelFunctorList.size() == m_PonderationList.size())
- {
- for (unsigned int i = 0; i<m_KernelFunctorList.size(); i++)
- {
- std::cout<<m_KernelFunctorList[i]->GetName()<<":"<<std::endl;
- std::cout<<"Associated ponderation:"<<m_PonderationList[i]<<std::endl;
- m_KernelFunctorList[i]->print_parameters();
- std::cout<<std::endl;
- }
- }
- else
- {
- itkGenericExceptionMacro(<<"ComposedKernelFunctor::print_param() : lists dimensions mismatch");
- }
-}
-
-
-int
-ComposedKernelFunctor::
-load_parameters(FILE ** pfile)
-{
- int NbParams(0);
- char keyword[81];
- char value[81];
-
- // Read functor name
- fscanf(*pfile,"%80s",keyword);
- this->SetName(std::string(keyword));
- // Read number of parameters
- fscanf(*pfile,"%d",&NbParams);
-
- for ( int cpt=0 ; cpt < NbParams ; cpt++)
- {
- fscanf(*pfile,"%80s",keyword);
- fscanf(*pfile,"%80s",value);
- this->SetValue<std::string>(keyword, value);
- }
-
- char tempChar[100];
- fscanf(*pfile, "%80s", tempChar);
- while( strcmp(tempChar,"Ponderation")==0 || strcmp(tempChar,"list:")==0 )
- {
- fscanf(*pfile, "%80s", tempChar);
- }
-
- unsigned int i = 0;
- while( strcmp(tempChar,"Kernels")!=0 )
- {
- m_PonderationList.push_back(::atof(tempChar));
- fscanf(*pfile, "%80s", tempChar);
- i++;
- }
- while( strcmp(tempChar,"Kernels")==0 || strcmp(tempChar,"list:")==0 || strcmp(tempChar,"Number")==0 || strcmp(tempChar,"of")==0 || strcmp(tempChar,"Kernels:")==0)
- {
- fscanf(*pfile, "%80s", tempChar);
- }
- int NbOfKernels = ::atoi(tempChar);
-
- for(unsigned int j=0; j<static_cast<unsigned int>(NbOfKernels); j++)
- {
- fscanf(*pfile, "%80s", tempChar);
- GenericKernelFunctorBase * gen;
- gen = new GenericKernelFunctorBase;
- gen->load_parameters(pfile);
- m_KernelFunctorList.push_back(gen);
- // Add the pointer to the "Have to Deleted" pointer list
- m_HaveToBeDeletedList.push_back(gen);
- }
-
- return 0;
-}
-
-
-int
-ComposedKernelFunctor::
-save_parameters(FILE ** pfile, const char * composed_kernel_parameters_keyword)const
-{
- MapConstIterator iter = this->GetMapParameters().begin();
- std::string line(composed_kernel_parameters_keyword);
- std::string strNbParams;
- ::otb::StringStream flux;
- flux << this->GetMapParameters().size();
- flux >> strNbParams;
- line = line + " " + this->GetName() + " " + strNbParams;
- while( iter != this->GetMapParameters().end() )
- {
- line = line + " " + iter->first + " " + iter->second;
- ++iter;
- }
- line = line + "\n" + "Ponderation list:\n";
-
- for (unsigned int i = 0; i<m_PonderationList.size(); i++)
- {
- std::string ponde;
- ::otb::StringStream flux;
- flux << m_PonderationList[i];
- flux >> ponde;
- line = line + " " + ponde;
- }
- line = line + "\n" + "Kernels list:\nNumber of Kernels: ";
- std::string nbOfKernels;
- ::otb::StringStream flux2;
- flux2 << m_KernelFunctorList.size();
- flux2 >> nbOfKernels;
-
- line = line + nbOfKernels + "\n";
- fprintf(*pfile,"%s", line.c_str());
- for (unsigned int i = 0; i<m_KernelFunctorList.size(); i++)
- {
- m_KernelFunctorList[i]->save_parameters(pfile, "generic_kernel_parameters");
- }
-
- return 0;
-}
-
-//} // namespace otb
diff --git a/Utilities/otbsvm/svm.h b/Utilities/otbsvm/svm.h
index e97eded539..168116212b 100644
--- a/Utilities/otbsvm/svm.h
+++ b/Utilities/otbsvm/svm.h
@@ -1,21 +1,22 @@
#ifndef _LIBSVM_H
#define _LIBSVM_H
-//OTB's modifications
-//namespace otb
-//{
-class GenericKernelFunctorBase;
-class ComposedKernelFunctor;
-//}
+#define LIBSVM_VERSION 300
+/*** Begin OTB modification ***/
#include <map>
#include <vector>
#include "otbMacro.h"
+class GenericKernelFunctorBase;
+class ComposedKernelFunctor;
+/*** End OTB modification ***/
#ifdef __cplusplus
extern "C" {
#endif
+extern int libsvm_version;
+
struct svm_node
{
int index;
@@ -30,210 +31,240 @@ struct svm_problem
};
enum { C_SVC, NU_SVC, ONE_CLASS, EPSILON_SVR, NU_SVR }; /* svm_type */
-//OTB's modifications
enum { LINEAR, POLY, RBF, SIGMOID, PRECOMPUTED, GENERIC, COMPOSED }; /* kernel_type */
struct svm_parameter
{
- int svm_type;
- int kernel_type;
- int degree; /* for poly */
- double gamma; /* for poly/rbf/sigmoid */
- double coef0; /* for poly/sigmoid */
- double const_coef; /* for otbSVMKernels.h. */
- double lin_coef; /* for otbSVMKernels.h. */
- char custom[500]; /* for user supplied kernel */
- //OTB's modifications : Use by the generic kernel
- /*otb::*/GenericKernelFunctorBase * kernel_generic;
- // Composed kernel
+ int svm_type;
+ int kernel_type;
+ int degree; /* for poly */
+ double gamma; /* for poly/rbf/sigmoid */
+ double coef0; /* for poly/sigmoid */
+
+ /*** Begin OTB modification ***/
+ double const_coef;
+ double lin_coef;
+ char custom[500];
+ GenericKernelFunctorBase * kernel_generic;
ComposedKernelFunctor * kernel_composed;
-
- /* these are for training only */
- double cache_size; /* in MB */
- double eps; /* stopping criteria */
- double C; /* for C_SVC, EPSILON_SVR and NU_SVR */
- int nr_weight; /* for C_SVC */
- int *weight_label; /* for C_SVC */
- double* weight; /* for C_SVC */
- double nu; /* for NU_SVC, ONE_CLASS, and NU_SVR */
- double p; /* for EPSILON_SVR */
- int shrinking; /* use the shrinking heuristics */
- int probability; /* do probability estimates */
+ /*** End OTB modification ***/
+
+ /* these are for training only */
+ double cache_size; /* in MB */
+ double eps; /* stopping criteria */
+ double C; /* for C_SVC, EPSILON_SVR and NU_SVR */
+ int nr_weight; /* for C_SVC */
+ int *weight_label; /* for C_SVC */
+ double* weight; /* for C_SVC */
+ double nu; /* for NU_SVC, ONE_CLASS, and NU_SVR */
+ double p; /* for EPSILON_SVR */
+ int shrinking; /* use the shrinking heuristics */
+ int probability; /* do probability estimates */
};
//
// svm_model
-//
+//
struct svm_model
{
- svm_parameter param; // parameter
- int nr_class; // number of classes, = 2 in regression/one class svm
- int l; // total #SV
- svm_node **SV; // SVs (SV[l])
- double **sv_coef; // coefficients for SVs in decision functions (sv_coef[k-1][l])
- double *rho; // constants in decision functions (rho[k*(k-1)/2])
- double *probA; // pariwise probability information
- double *probB;
-
- // for classification only
-
- int *label; // label of each class (label[k])
- int *nSV; // number of SVs for each class (nSV[k])
- // nSV[0] + nSV[1] + ... + nSV[k-1] = l
- // XXX
- int free_sv; // 1 if svm_model is created by svm_load_model
+ struct svm_parameter param; /* parameter */
+ int nr_class; /* number of classes, = 2 in regression/one class svm */
+ int l; /* total #SV */
+ struct svm_node **SV; /* SVs (SV[l]) */
+ double **sv_coef; /* coefficients for SVs in decision functions (sv_coef[k-1][l]) */
+ double *rho; /* constants in decision functions (rho[k*(k-1)/2]) */
+ double *probA; /* pariwise probability information */
+ double *probB;
+
+ /* for classification only */
+
+ int *label; /* label of each class (label[k]) */
+ int *nSV; /* number of SVs for each class (nSV[k]) */
+ /* nSV[0] + nSV[1] + ... + nSV[k-1] = l */
+ /* XXX */
+ int free_sv; /* 1 if svm_model is created by svm_load_model*/
+ /* 0 if svm_model is created by svm_train */
+
+ /*** Begin OTB modification ***/
bool delete_composed; // to know if the composed functor was set using load method
- // 0 if svm_model is created by svm_train
+ /*** End OTB modification ***/
};
struct svm_model *svm_train(const struct svm_problem *prob, const struct svm_parameter *param);
void svm_cross_validation(const struct svm_problem *prob, const struct svm_parameter *param, int nr_fold, double *target);
int svm_save_model(const char *model_file_name, const struct svm_model *model);
+/*** Begin OTB modification ***/
+struct svm_model *svm_load_model(const char *model_file_name, GenericKernelFunctorBase* generic_kernel_functor = NULL);
+struct svm_model *svm_copy_model( const svm_model *model );
+/*** End OTB modification ***/
int svm_get_svm_type(const struct svm_model *model);
int svm_get_nr_class(const struct svm_model *model);
void svm_get_labels(const struct svm_model *model, int *label);
double svm_get_svr_probability(const struct svm_model *model);
-void svm_predict_values(const struct svm_model *model, const struct svm_node *x, double* dec_values);
+double svm_predict_values(const struct svm_model *model, const struct svm_node *x, double* dec_values);
double svm_predict(const struct svm_model *model, const struct svm_node *x);
double svm_predict_probability(const struct svm_model *model, const struct svm_node *x, double* prob_estimates);
-void svm_destroy_model(struct svm_model *model);
+void svm_free_model_content(struct svm_model *model_ptr);
+void svm_free_and_destroy_model(struct svm_model **model_ptr_ptr);
void svm_destroy_param(struct svm_parameter *param);
const char *svm_check_parameter(const struct svm_problem *prob, const struct svm_parameter *param);
int svm_check_probability_model(const struct svm_model *model);
-//OTB's modifications
-struct svm_model *svm_load_model(const char *model_file_name, /*otb::*/GenericKernelFunctorBase * generic_kernel_functor = NULL);
-struct svm_model *svm_copy_model( const svm_model *model );
+void svm_set_print_string_function(void (*print_func)(const char *));
+
+// deprecated
+// this function will be removed in future release
+void svm_destroy_model(struct svm_model *model_ptr);
#ifdef __cplusplus
}
#endif
+/*** Begin OTB modification ***/
-/*
-namespace otb
-{
-*/
-//OTB's modifications
/** \class GenericKernelFunctorBase
* \brief Undocumented
*/
class GenericKernelFunctorBase
{
public:
- GenericKernelFunctorBase() : m_Name("FunctorName") {};
+ GenericKernelFunctorBase() :
+ m_Name("FunctorName")
+ {
+ }
/** Recopy constructor */
- GenericKernelFunctorBase( const GenericKernelFunctorBase& copy);
+ GenericKernelFunctorBase(const GenericKernelFunctorBase& copy);
virtual GenericKernelFunctorBase& operator=(const GenericKernelFunctorBase& copy);
- virtual ~GenericKernelFunctorBase() {};
+ virtual ~GenericKernelFunctorBase()
+ {
+ }
- typedef GenericKernelFunctorBase Superclass;
- typedef std::map<std::string,std::string> MapType;
- typedef MapType::iterator MapIterator;
- typedef MapType::const_iterator MapConstIterator;
+ typedef GenericKernelFunctorBase Superclass;
+ typedef std::map<std::string, std::string> MapType;
+ typedef MapType::iterator MapIterator;
+ typedef MapType::const_iterator MapConstIterator;
template<class T>
T GetValue(const char *option) const
{
- std::string Value = m_MapParameters.find(std::string(option))->second;
- T lValeur;
- ::otb::StringStream flux;
- flux << Value;
- flux >> lValeur;
- return lValeur;
+ std::string Value = m_MapParameters.find(std::string(option))->second;
+ T lValeur;
+ ::otb::StringStream flux;
+ flux << Value;
+ flux >> lValeur;
+ return lValeur;
}
+
template<class T>
void SetValue(const char *option, const T & value)
{
- std::string lValeur;
- ::otb::StringStream flux;
- flux << value;
- flux >> lValeur;
- m_MapParameters[std::string(option)] = lValeur;
+ std::string lValeur;
+ ::otb::StringStream flux;
+ flux << value;
+ flux >> lValeur;
+ m_MapParameters[std::string(option)] = lValeur;
}
-
- virtual double operator()(const svm_node * /*x*/, const svm_node * /*y*/, const svm_parameter& /*param*/)const
- {
- itkGenericExceptionMacro(<<"Kernel functor not definied (Null)");
- return static_cast<double>(0.);
- }
+ virtual double operator()(const svm_node * /*x*/, const svm_node * /*y*/, const svm_parameter& /*param*/) const
+ {
+ itkGenericExceptionMacro(<<"Kernel functor not definied (Null)");
+ return static_cast<double> (0.);
+ }
/** Used for Taylor classification*/
// degree is the development degree
// index is the current value
// isAtEnd to indicate that it's the last possible derivation
// baseValue is the constant of the formula
- virtual double derivative(const svm_node * /*x*/, const svm_node * /*y*/, const svm_parameter& /*param*/, int /*degree*/, int /*index*/, bool /*isAtEnd*/, double /*constValue*/)const
- {
- itkGenericExceptionMacro(<<"derivative method not definied (Null)");
- return 0.;
- }
+ virtual double derivative(const svm_node * /*x*/, const svm_node * /*y*/, const svm_parameter& /*param*/,
+ int /*degree*/, int /*index*/, bool /*isAtEnd*/, double /*constValue*/) const
+ {
+ itkGenericExceptionMacro(<<"derivative method not definied (Null)");
+ return 0.;
+ }
virtual int load_parameters(FILE ** pfile);
- virtual int save_parameters(FILE ** pfile, const char * generic_kernel_parameters_keyword)const;
+ virtual int save_parameters(FILE ** pfile, const char * generic_kernel_parameters_keyword) const;
- virtual void print_parameters(void)const;
+ virtual void print_parameters(void) const;
- virtual void Update(void){ }
+ virtual double dot(const svm_node *px, const svm_node *py) const;
- virtual double dot(const svm_node *px, const svm_node *py)const;
+ virtual svm_node * sub(const svm_node *px, const svm_node *py) const;
- virtual svm_node * sub(const svm_node *px, const svm_node *py)const;
+ virtual svm_node * add(const svm_node *px, const svm_node *py) const;
- virtual svm_node * add(const svm_node *px, const svm_node *py)const;
+ virtual void SetName(std::string name)
+ {
+ m_Name = name;
+ }
+ virtual std::string GetName(void)
+ {
+ return m_Name;
+ }
+ virtual const std::string GetName(void) const
+ {
+ return m_Name;
+ }
- virtual void SetName(std::string name ) { m_Name = name;}
- virtual std::string GetName(void) { return m_Name;}
- virtual const std::string GetName(void) const { return m_Name;}
+ virtual void SetMapParameters(const MapType & map)
+ {
+ m_MapParameters = map;
+ }
+
+ virtual const MapType & GetMapParameters() const
+ {
+ return m_MapParameters;
+ }
- virtual void SetMapParameters(const MapType & map){ m_MapParameters = map; };
- virtual const MapType & GetMapParameters()const { return m_MapParameters; };
- virtual MapType GetMapParameters(){ return m_MapParameters; };
+ virtual MapType GetMapParameters()
+ {
+ return m_MapParameters;
+ }
private:
- /** Kernel functor parameters */
- MapType m_MapParameters;
+ /** Kernel functor parameters */
+ MapType m_MapParameters;
- /** Functor label name (without space) */
- std::string m_Name;
+ /** Functor label name (without space) */
+ std::string m_Name;
};
-
/** \class ComposedKernelFunctor
* \brief Undocumented
*/
-class ComposedKernelFunctor : public GenericKernelFunctorBase
+class ComposedKernelFunctor: public GenericKernelFunctorBase
{
public:
ComposedKernelFunctor()
- {
- this->SetName("ComposedFunctorName");
- this->SetValue<bool>("MultiplyKernelFunctor", false);
- };
+ {
+ this->SetName("ComposedFunctorName");
+ this->SetValue<bool> ("MultiplyKernelFunctor", false);
+ }
+
virtual ~ComposedKernelFunctor()
- {
- for(unsigned int i=0; i<m_HaveToBeDeletedList.size(); i++)
- {
- for(unsigned int j=0; j<m_KernelFunctorList.size(); j++)
- {
- if(m_KernelFunctorList[j] == m_HaveToBeDeletedList[i])
- {
- delete m_KernelFunctorList[j];
- m_HaveToBeDeletedList[i] = NULL;
- }
- }
- }
- };
+ {
+ for(unsigned int i=0; i < m_HaveToBeDeletedList.size(); i++)
+ {
+ for(unsigned int j=0; j<m_KernelFunctorList.size(); j++)
+ {
+ if(m_KernelFunctorList[j] == m_HaveToBeDeletedList[i])
+ {
+ delete m_KernelFunctorList[j];
+ m_HaveToBeDeletedList[i] = NULL;
+ }
+ }
+ }
+ }
/** Recopy constructor */
ComposedKernelFunctor( const ComposedKernelFunctor& copy );
@@ -241,28 +272,28 @@ public:
typedef std::vector<GenericKernelFunctorBase *> KernelListType;
- virtual double operator()(const svm_node *x, const svm_node *y, const svm_parameter& param)const // = 0
+ virtual double operator()(const svm_node *x, const svm_node *y, const svm_parameter& param) const
{
- double out = 0.;
- if (m_KernelFunctorList.size() != 0 && m_PonderationList.size() != 0 && m_KernelFunctorList.size() == m_PonderationList.size())
- {
- for (unsigned int i = 0; i<m_KernelFunctorList.size(); i++)
- {
- if ((this->GetValue<bool>("MultiplyKernelFunctor")) == false)
- {
- out += m_PonderationList[i]*(*m_KernelFunctorList[i])(x, y, param);
- }
- else
- {
- out *= (*m_KernelFunctorList[i])(x, y, param);
- }
- }
- }
- else
- {
- itkGenericExceptionMacro(<<"ComposedKernelFunctor::operator() : lists dimensions mismatch");
- }
- return out;
+ double out = 0.;
+ if (m_KernelFunctorList.size() != 0 && m_PonderationList.size() != 0 && m_KernelFunctorList.size() == m_PonderationList.size())
+ {
+ for (unsigned int i = 0; i<m_KernelFunctorList.size(); i++)
+ {
+ if ((this->GetValue<bool>("MultiplyKernelFunctor")) == false)
+ {
+ out += m_PonderationList[i]*(*m_KernelFunctorList[i])(x, y, param);
+ }
+ else
+ {
+ out *= (*m_KernelFunctorList[i])(x, y, param);
+ }
+ }
+ }
+ else
+ {
+ itkGenericExceptionMacro(<<"ComposedKernelFunctor::operator() : lists dimensions mismatch");
+ }
+ return out;
}
/** Used for Taylor classification*/
@@ -271,28 +302,28 @@ public:
// isAtEnd to indicate that it's the last possible derivation
// baseValue is the constant of the formula
virtual double derivative(const svm_node *x, const svm_node *y, const svm_parameter& param, int degree, int index, bool isAtEnd, double constValue)const
- {
- double out = 0.;
- if (m_KernelFunctorList.size() != 0 && m_PonderationList.size() != 0 && m_KernelFunctorList.size() == m_PonderationList.size())
- {
- for (unsigned int i = 0; i<m_KernelFunctorList.size(); i++)
- {
- if ((this->GetValue<bool>("MultiplyKernelFunctor")) == false)
- {
- out += m_PonderationList[i]*(m_KernelFunctorList[i]->derivative(x, y, param, degree, index, isAtEnd, constValue));
- }
- else
- {
- itkGenericExceptionMacro(<<"derivative method not definied (Null)");
- }
- }
- }
- else
- {
- itkGenericExceptionMacro(<<"ComposedKernelFunctor::operator() : lists dimensions mismatch");
- }
- return out;
- }
+ {
+ double out = 0.;
+ if (m_KernelFunctorList.size() != 0 && m_PonderationList.size() != 0 && m_KernelFunctorList.size() == m_PonderationList.size())
+ {
+ for (unsigned int i = 0; i<m_KernelFunctorList.size(); i++)
+ {
+ if ((this->GetValue<bool>("MultiplyKernelFunctor")) == false)
+ {
+ out += m_PonderationList[i]*(m_KernelFunctorList[i]->derivative(x, y, param, degree, index, isAtEnd, constValue));
+ }
+ else
+ {
+ itkGenericExceptionMacro(<<"derivative method not definied (Null)");
+ }
+ }
+ }
+ else
+ {
+ itkGenericExceptionMacro(<<"ComposedKernelFunctor::operator() : lists dimensions mismatch");
+ }
+ return out;
+ }
virtual int load_parameters(FILE ** pfile);
@@ -300,50 +331,56 @@ public:
virtual void print_parameters(void)const;
- //virtual void Update(void){};
+ /** Set/Get the SVM Model vector for the composed kernel */
+ KernelListType GetKernelFunctorList()
+ {return m_KernelFunctorList;};
+ void SetKernelFunctorList(KernelListType kernelFunctorList)
+ {m_KernelFunctorList = kernelFunctorList;};
- /** Set/Get the SVM Model vector for the composed kernel */
- KernelListType GetKernelFunctorList(){ return m_KernelFunctorList; };
- void SetKernelFunctorList(KernelListType kernelFunctorList){ m_KernelFunctorList = kernelFunctorList; };
// Add 1 element to the end of the list
- void AddKernelFunctorModelToKernelList(GenericKernelFunctorBase * kernelfunctor){ m_KernelFunctorList.push_back(kernelfunctor); };
+ void AddKernelFunctorModelToKernelList(GenericKernelFunctorBase * kernelfunctor)
+ {m_KernelFunctorList.push_back(kernelfunctor);};
/** Generic kernel functors that have to be deleted. */
- KernelListType GetHaveToBeDeletedList(){ return m_HaveToBeDeletedList; };
- void SetHaveToBeDeletedList(KernelListType kernelFunctorList){ m_HaveToBeDeletedList = kernelFunctorList; };
+ KernelListType GetHaveToBeDeletedList()
+ {return m_HaveToBeDeletedList;};
+ void SetHaveToBeDeletedList(KernelListType kernelFunctorList)
+ {m_HaveToBeDeletedList = kernelFunctorList;};
// Add 1 element to the end of the list
- void AddKernelFunctorModelToDeleteKernelList(GenericKernelFunctorBase * kernelfunctor){ m_HaveToBeDeletedList.push_back(kernelfunctor); };
+ void AddKernelFunctorModelToDeleteKernelList(GenericKernelFunctorBase * kernelfunctor)
+ {m_HaveToBeDeletedList.push_back(kernelfunctor);};
/** Set/Get the ponderation list to apply to each svm_model of the composed kernel */
- std::vector<double> GetPonderationList(){ return m_PonderationList; };
- void SetPonderationModelList(const std::vector<double> & list){ m_PonderationList = list; };
+ std::vector<double> GetPonderationList()
+ {return m_PonderationList;};
+ void SetPonderationModelList(const std::vector<double> & list)
+ {m_PonderationList = list;};
// Add 1 element to the end of the list
- void AddPonderationToPonderationList(const double & pond){ m_PonderationList.push_back(pond); };
+ void AddPonderationToPonderationList(const double & pond)
+ {m_PonderationList.push_back(pond);};
- /** Set/Get the boolean to know which operation has to be done with the kernel functors. */
- void SetMultiplyKernelFunctor( bool val ){ this->SetValue<bool>("MultiplyKernelFunctor", val); };
- bool GetMultiplyKernelFunctor(){ return (this->GetValue<bool>("MultiplyKernelFunctor")); };
+ /** Set/Get the boolean to know which operation has to be done with the kernel functors. */
+ void SetMultiplyKernelFunctor( bool val )
+ {this->SetValue<bool>("MultiplyKernelFunctor", val);};
+ bool GetMultiplyKernelFunctor()
+ {return (this->GetValue<bool>("MultiplyKernelFunctor"));};
private:
- typedef GenericKernelFunctorBase::MapType MapType;
- typedef GenericKernelFunctorBase::MapIterator MapIterator;
- typedef GenericKernelFunctorBase::MapConstIterator MapConstIterator;
-
+ typedef GenericKernelFunctorBase::MapType MapType;
+ typedef GenericKernelFunctorBase::MapIterator MapIterator;
+ typedef GenericKernelFunctorBase::MapConstIterator MapConstIterator;
- /** Generic kernel functors that composed kernel */
- KernelListType m_KernelFunctorList;
+ /** Generic kernel functors that composed kernel */
+ KernelListType m_KernelFunctorList;
/** Generic kernel functors that have to be deleted.
- * This list was made for the load_parameters methods where you set new functors using new.
- * But, in other cases, functor can be added with reference. Thus, we need to know which ones have to be deleted. */
- KernelListType m_HaveToBeDeletedList;
- /** Ponderation list to apply to each svm_model of the composed kernel*/
- std::vector<double> m_PonderationList;
-
+ * This list was made for the load_parameters methods where you set new functors using new.
+ * But, in other cases, functor can be added with reference. Thus, we need to know which ones have to be deleted. */
+ KernelListType m_HaveToBeDeletedList;
+ /** Ponderation list to apply to each svm_model of the composed kernel*/
+ std::vector<double> m_PonderationList;
};
-
-
-//} // namespace otb
+/*** End OTB modification ***/
#endif /* _LIBSVM_H */
--
GitLab