[pktools] 18/375: added pkclassify_svm.cc and pkdsm2shadow.cc
Bas Couwenberg
sebastic at xs4all.nl
Wed Dec 3 21:53:54 UTC 2014
This is an automated email from the git hooks/post-receive script.
sebastic-guest pushed a commit to branch upstream-master
in repository pktools.
commit bd1c222b76545fdb5dc310736f1f4fb06b49de68
Author: Pieter Kempeneers <kempenep at gmail.com>
Date: Wed Oct 31 18:07:25 2012 +0100
added pkclassify_svm.cc and pkdsm2shadow.cc
---
README | 8 +-
src/algorithms/SVM_COPYRIGHT | 31 +
src/algorithms/svm.cpp | 3114 +++++++++++++++++++++++++++++++++++++
src/algorithms/svm.h | 101 ++
src/apps/Makefile.am | 4 +-
src/apps/Makefile.in | 48 +-
src/apps/pkclassify_svm.cc | 1037 ++++++++++++
src/apps/pkdsm2shadow.cc | 148 ++
src/apps/pkextract.cc | 4 +-
src/imageclasses/ImgWriterGdal.cc | 2 +
10 files changed, 4466 insertions(+), 31 deletions(-)
diff --git a/README b/README
index 12bcc9c..aa8b4ca 100644
--- a/README
+++ b/README
@@ -14,7 +14,7 @@ To install the programs in pktools, refer to the file INSTALL
Change history
-------------
-June 25 2012, first public release of the code
-September 04 2012, introduced --enable-fann and --enable-las in configuration
-September 13 2012, support spectral filtering (z-dimension) in pkfilter using tapz option
-October 20 2012, support for SVM classifier
+version 1.0 June 25 2012, first public release of the code
+version 2.1 September 04 2012, introduced --enable-fann and --enable-las in configuration
+version 2.2 September 13 2012, support spectral filtering (z-dimension) in pkfilter using tapz option
+version 2.3 October 20 2012, support for SVM classifier
diff --git a/src/algorithms/SVM_COPYRIGHT b/src/algorithms/SVM_COPYRIGHT
new file mode 100644
index 0000000..9280fc0
--- /dev/null
+++ b/src/algorithms/SVM_COPYRIGHT
@@ -0,0 +1,31 @@
+
+Copyright (c) 2000-2012 Chih-Chung Chang and Chih-Jen Lin
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+1. Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright
+notice, this list of conditions and the following disclaimer in the
+documentation and/or other materials provided with the distribution.
+
+3. Neither name of copyright holders nor the names of its contributors
+may be used to endorse or promote products derived from this software
+without specific prior written permission.
+
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR
+CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/src/algorithms/svm.cpp b/src/algorithms/svm.cpp
new file mode 100644
index 0000000..5e9fd28
--- /dev/null
+++ b/src/algorithms/svm.cpp
@@ -0,0 +1,3114 @@
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <ctype.h>
+#include <float.h>
+#include <string.h>
+#include <stdarg.h>
+#include <limits.h>
+#include <locale.h>
+//test
+// #include <iostream>
+#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:
+
+ double (Kernel::*kernel_function)(int i, int j) const;
+
+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
+ {
+ return dot(x[i],x[j]);
+ }
+ double kernel_poly(int i, int j) const
+ {
+ return powi(gamma*dot(x[i],x[j])+coef0,degree);
+ }
+ double kernel_rbf(int i, int j) const
+ {
+ return exp(-gamma*(x_square[i]+x_square[j]-2*dot(x[i],x[j])));
+ }
+ double kernel_sigmoid(int i, int j) const
+ {
+ return tanh(gamma*dot(x[i],x[j])+coef0);
+ }
+ double kernel_precomputed(int i, int j) const
+ {
+ return x[i][(int)(x[j][0].value)].value;
+ }
+};
+
+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)
+{
+ 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;
+ }
+
+ 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;
+ 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 max_iter = max(10000000, l>INT_MAX/100 ? INT_MAX : 100*l);
+ int counter = min(l,1000)+1;
+
+ while(iter < max_iter)
+ {
+ // 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];
+ }
+ }
+ }
+
+ if(iter >= max_iter)
+ {
+ if(active_size < l)
+ {
+ // reconstruct the whole gradient to calculate objective value
+ reconstruct_gradient();
+ active_size = l;
+ info("*");
+ }
+ info("\nWARNING: reaching max number of iterations");
+ }
+
+ // 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++)
+ QD[i] = (this->*kernel_function)(i,i);
+ }
+
+ 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++)
+ data[j] = (Qfloat)(y[i]*y[j]*(this->*kernel_function)(i,j));
+ }
+ 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++)
+ QD[i] = (this->*kernel_function)(i,i);
+ }
+
+ 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++)
+ data[j] = (Qfloat)(this->*kernel_function)(i,j);
+ }
+ 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;
+ QD[k] = (this->*kernel_function)(k,k);
+ 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++)
+ data[j] = (Qfloat)(this->*kernel_function)(real_i,j);
+ }
+
+ // 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;
+ // 1-p used later; avoid catastrophic cancellation
+ 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);
+ if(nr_class == 1)
+ info("WARNING: training data in only one class. See README for details.\n");
+
+ 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)
+{
+ int i;
+ 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(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 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[]=
+{
+ "linear","polynomial","rbf","sigmoid","precomputed",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;
+
+ char *old_locale = strdup(setlocale(LC_ALL, NULL));
+ setlocale(LC_ALL, "C");
+
+ 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]);
+
+ 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");
+ }
+
+ setlocale(LC_ALL, old_locale);
+ free(old_locale);
+
+ 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)
+{
+ FILE *fp = fopen(model_file_name,"rb");
+ if(fp==NULL) return NULL;
+
+ char *old_locale = strdup(setlocale(LC_ALL, NULL));
+ setlocale(LC_ALL, "C");
+
+ // 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");
+
+ setlocale(LC_ALL, old_locale);
+ free(old_locale);
+ 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");
+
+ setlocale(LC_ALL, old_locale);
+ free(old_locale);
+ 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;
+ }
+ else
+ {
+ fprintf(stderr,"unknown text in model file: [%s]\n",cmd);
+
+ setlocale(LC_ALL, old_locale);
+ free(old_locale);
+ 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);
+
+ setlocale(LC_ALL, old_locale);
+ free(old_locale);
+
+ if (ferror(fp) != 0 || fclose(fp) != 0)
+ return NULL;
+
+ model->free_sv = 1; // XXX
+ return model;
+}
+
+void svm_free_model_content(svm_model* model_ptr)
+{
+ if(model_ptr->free_sv && model_ptr->l > 0 && model_ptr->SV != NULL)
+ free((void *)(model_ptr->SV[0]));
+ if(model_ptr->sv_coef)
+ {
+ for(int i=0;i<model_ptr->nr_class-1;i++)
+ free(model_ptr->sv_coef[i]);
+ }
+
+ free(model_ptr->SV);
+ model_ptr->SV = NULL;
+
+ free(model_ptr->sv_coef);
+ model_ptr->sv_coef = NULL;
+
+ free(model_ptr->rho);
+ model_ptr->rho = NULL;
+
+ free(model_ptr->label);
+ model_ptr->label= NULL;
+
+ free(model_ptr->probA);
+ model_ptr->probA = NULL;
+
+ free(model_ptr->probB);
+ model_ptr->probB= NULL;
+
+ free(model_ptr->nSV);
+ model_ptr->nSV = NULL;
+}
+
+void svm_free_and_destroy_model(svm_model** model_ptr_ptr)
+{
+ if(model_ptr_ptr != NULL && *model_ptr_ptr != NULL)
+ {
+ svm_free_model_content(*model_ptr_ptr);
+ free(*model_ptr_ptr);
+ *model_ptr_ptr = NULL;
+ }
+}
+
+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)
+ 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;
+}
diff --git a/src/algorithms/svm.h b/src/algorithms/svm.h
new file mode 100644
index 0000000..2f60a57
--- /dev/null
+++ b/src/algorithms/svm.h
@@ -0,0 +1,101 @@
+#ifndef _LIBSVM_H
+#define _LIBSVM_H
+
+#define LIBSVM_VERSION 312
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+extern int libsvm_version;
+
+struct svm_node
+{
+ int index;
+ double value;
+};
+
+struct svm_problem
+{
+ int l;
+ double *y;
+ struct svm_node **x;
+};
+
+enum { C_SVC, NU_SVC, ONE_CLASS, EPSILON_SVR, NU_SVR }; /* svm_type */
+enum { LINEAR, POLY, RBF, SIGMOID, PRECOMPUTED }; /* 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 */
+
+ /* 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
+{
+ 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 */
+};
+
+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);
+struct svm_model *svm_load_model(const char *model_file_name);
+
+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);
+
+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_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);
+
+void svm_set_print_string_function(void (*print_func)(const char *));
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* _LIBSVM_H */
diff --git a/src/apps/Makefile.am b/src/apps/Makefile.am
index 3734a79..d8b1238 100644
--- a/src/apps/Makefile.am
+++ b/src/apps/Makefile.am
@@ -6,7 +6,7 @@ LDADD = $(GDAL_LDFLAGS) $(top_builddir)/src/algorithms/libalgorithms.a $(top_bui
###############################################################################
# the program to build (the names of the final binaries)
-bin_PROGRAMS = pkinfo pkcrop pkreclass pkgetmask pksetmask pkcreatect pkdumpimg pkdumpogr pksieve pkstat pkstatogr pkegcs pkextract pkfillnodata pkfilter pkveg2shadow pkmosaic pkndvi pkpolygonize pkascii2img pkdiff pkclassify_svm
+bin_PROGRAMS = pkinfo pkcrop pkreclass pkgetmask pksetmask pkcreatect pkdumpimg pkdumpogr pksieve pkstat pkstatogr pkegcs pkextract pkfillnodata pkfilter pkdsm2shadow pkmosaic pkndvi pkpolygonize pkascii2img pkdiff pkclassify_svm
if USE_FANN
bin_PROGRAMS += pkclassify_nn
pkclassify_nn_SOURCES = $(top_srcdir)/src/algorithms/myfann_cpp.h pkclassify_nn.h pkclassify_nn.cc
@@ -38,7 +38,7 @@ pkextract_SOURCES = pkextract.cc
pkfillnodata_SOURCES = pkfillnodata.cc
pkfilter_SOURCES = pkfilter.cc
#pkfilter_LDADD = -lgdal $(AM_LDFLAGS) -lgdal
-pkveg2shadow_SOURCES = pkveg2shadow.cc
+pkdsm2shadow_SOURCES = pkdsm2shadow.cc
pkmosaic_SOURCES = pkmosaic.cc
pkndvi_SOURCES = pkndvi.cc
pkpolygonize_SOURCES = pkpolygonize.cc
diff --git a/src/apps/Makefile.in b/src/apps/Makefile.in
index 4a8dbbe..f83de84 100644
--- a/src/apps/Makefile.in
+++ b/src/apps/Makefile.in
@@ -37,7 +37,7 @@ bin_PROGRAMS = pkinfo$(EXEEXT) pkcrop$(EXEEXT) pkreclass$(EXEEXT) \
pkdumpimg$(EXEEXT) pkdumpogr$(EXEEXT) pksieve$(EXEEXT) \
pkstat$(EXEEXT) pkstatogr$(EXEEXT) pkegcs$(EXEEXT) \
pkextract$(EXEEXT) pkfillnodata$(EXEEXT) pkfilter$(EXEEXT) \
- pkveg2shadow$(EXEEXT) pkmosaic$(EXEEXT) pkndvi$(EXEEXT) \
+ pkdsm2shadow$(EXEEXT) pkmosaic$(EXEEXT) pkndvi$(EXEEXT) \
pkpolygonize$(EXEEXT) pkascii2img$(EXEEXT) pkdiff$(EXEEXT) \
pkclassify_svm$(EXEEXT) $(am__EXEEXT_1) $(am__EXEEXT_2)
@USE_FANN_TRUE at am__append_1 = pkclassify_nn
@@ -101,6 +101,12 @@ pkdiff_LDADD = $(LDADD)
pkdiff_DEPENDENCIES = $(am__DEPENDENCIES_1) \
$(top_builddir)/src/algorithms/libalgorithms.a \
$(top_builddir)/src/imageclasses/libimageClasses.a
+am_pkdsm2shadow_OBJECTS = pkdsm2shadow.$(OBJEXT)
+pkdsm2shadow_OBJECTS = $(am_pkdsm2shadow_OBJECTS)
+pkdsm2shadow_LDADD = $(LDADD)
+pkdsm2shadow_DEPENDENCIES = $(am__DEPENDENCIES_1) \
+ $(top_builddir)/src/algorithms/libalgorithms.a \
+ $(top_builddir)/src/imageclasses/libimageClasses.a
am_pkdumpimg_OBJECTS = pkdumpimg.$(OBJEXT)
pkdumpimg_OBJECTS = $(am_pkdumpimg_OBJECTS)
pkdumpimg_LDADD = $(LDADD)
@@ -195,12 +201,6 @@ pkstatogr_LDADD = $(LDADD)
pkstatogr_DEPENDENCIES = $(am__DEPENDENCIES_1) \
$(top_builddir)/src/algorithms/libalgorithms.a \
$(top_builddir)/src/imageclasses/libimageClasses.a
-am_pkveg2shadow_OBJECTS = pkveg2shadow.$(OBJEXT)
-pkveg2shadow_OBJECTS = $(am_pkveg2shadow_OBJECTS)
-pkveg2shadow_LDADD = $(LDADD)
-pkveg2shadow_DEPENDENCIES = $(am__DEPENDENCIES_1) \
- $(top_builddir)/src/algorithms/libalgorithms.a \
- $(top_builddir)/src/imageclasses/libimageClasses.a
DEFAULT_INCLUDES = -I. at am__isrc@ -I$(top_builddir)
depcomp = $(SHELL) $(top_srcdir)/depcomp
am__depfiles_maybe = depfiles
@@ -216,23 +216,25 @@ CCLD = $(CC)
LINK = $(CCLD) $(AM_CFLAGS) $(CFLAGS) $(AM_LDFLAGS) $(LDFLAGS) -o $@
SOURCES = $(pkascii2img_SOURCES) $(pkclassify_nn_SOURCES) \
$(pkclassify_svm_SOURCES) $(pkcreatect_SOURCES) \
- $(pkcrop_SOURCES) $(pkdiff_SOURCES) $(pkdumpimg_SOURCES) \
- $(pkdumpogr_SOURCES) $(pkegcs_SOURCES) $(pkextract_SOURCES) \
- $(pkfillnodata_SOURCES) $(pkfilter_SOURCES) \
- $(pkgetmask_SOURCES) $(pkinfo_SOURCES) $(pklas2img_SOURCES) \
- $(pkmosaic_SOURCES) $(pkndvi_SOURCES) $(pkpolygonize_SOURCES) \
- $(pkreclass_SOURCES) $(pksetmask_SOURCES) $(pksieve_SOURCES) \
- $(pkstat_SOURCES) $(pkstatogr_SOURCES) $(pkveg2shadow_SOURCES)
-DIST_SOURCES = $(pkascii2img_SOURCES) \
- $(am__pkclassify_nn_SOURCES_DIST) $(pkclassify_svm_SOURCES) \
- $(pkcreatect_SOURCES) $(pkcrop_SOURCES) $(pkdiff_SOURCES) \
+ $(pkcrop_SOURCES) $(pkdiff_SOURCES) $(pkdsm2shadow_SOURCES) \
$(pkdumpimg_SOURCES) $(pkdumpogr_SOURCES) $(pkegcs_SOURCES) \
$(pkextract_SOURCES) $(pkfillnodata_SOURCES) \
$(pkfilter_SOURCES) $(pkgetmask_SOURCES) $(pkinfo_SOURCES) \
+ $(pklas2img_SOURCES) $(pkmosaic_SOURCES) $(pkndvi_SOURCES) \
+ $(pkpolygonize_SOURCES) $(pkreclass_SOURCES) \
+ $(pksetmask_SOURCES) $(pksieve_SOURCES) $(pkstat_SOURCES) \
+ $(pkstatogr_SOURCES)
+DIST_SOURCES = $(pkascii2img_SOURCES) \
+ $(am__pkclassify_nn_SOURCES_DIST) $(pkclassify_svm_SOURCES) \
+ $(pkcreatect_SOURCES) $(pkcrop_SOURCES) $(pkdiff_SOURCES) \
+ $(pkdsm2shadow_SOURCES) $(pkdumpimg_SOURCES) \
+ $(pkdumpogr_SOURCES) $(pkegcs_SOURCES) $(pkextract_SOURCES) \
+ $(pkfillnodata_SOURCES) $(pkfilter_SOURCES) \
+ $(pkgetmask_SOURCES) $(pkinfo_SOURCES) \
$(am__pklas2img_SOURCES_DIST) $(pkmosaic_SOURCES) \
$(pkndvi_SOURCES) $(pkpolygonize_SOURCES) $(pkreclass_SOURCES) \
$(pksetmask_SOURCES) $(pksieve_SOURCES) $(pkstat_SOURCES) \
- $(pkstatogr_SOURCES) $(pkveg2shadow_SOURCES)
+ $(pkstatogr_SOURCES)
ETAGS = etags
CTAGS = ctags
DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST)
@@ -371,7 +373,7 @@ pkextract_SOURCES = pkextract.cc
pkfillnodata_SOURCES = pkfillnodata.cc
pkfilter_SOURCES = pkfilter.cc
#pkfilter_LDADD = -lgdal $(AM_LDFLAGS) -lgdal
-pkveg2shadow_SOURCES = pkveg2shadow.cc
+pkdsm2shadow_SOURCES = pkdsm2shadow.cc
pkmosaic_SOURCES = pkmosaic.cc
pkndvi_SOURCES = pkndvi.cc
pkpolygonize_SOURCES = pkpolygonize.cc
@@ -467,6 +469,9 @@ pkcrop$(EXEEXT): $(pkcrop_OBJECTS) $(pkcrop_DEPENDENCIES)
pkdiff$(EXEEXT): $(pkdiff_OBJECTS) $(pkdiff_DEPENDENCIES)
@rm -f pkdiff$(EXEEXT)
$(CXXLINK) $(pkdiff_OBJECTS) $(pkdiff_LDADD) $(LIBS)
+pkdsm2shadow$(EXEEXT): $(pkdsm2shadow_OBJECTS) $(pkdsm2shadow_DEPENDENCIES)
+ @rm -f pkdsm2shadow$(EXEEXT)
+ $(CXXLINK) $(pkdsm2shadow_OBJECTS) $(pkdsm2shadow_LDADD) $(LIBS)
pkdumpimg$(EXEEXT): $(pkdumpimg_OBJECTS) $(pkdumpimg_DEPENDENCIES)
@rm -f pkdumpimg$(EXEEXT)
$(CXXLINK) $(pkdumpimg_OBJECTS) $(pkdumpimg_LDADD) $(LIBS)
@@ -518,9 +523,6 @@ pkstat$(EXEEXT): $(pkstat_OBJECTS) $(pkstat_DEPENDENCIES)
pkstatogr$(EXEEXT): $(pkstatogr_OBJECTS) $(pkstatogr_DEPENDENCIES)
@rm -f pkstatogr$(EXEEXT)
$(CXXLINK) $(pkstatogr_OBJECTS) $(pkstatogr_LDADD) $(LIBS)
-pkveg2shadow$(EXEEXT): $(pkveg2shadow_OBJECTS) $(pkveg2shadow_DEPENDENCIES)
- @rm -f pkveg2shadow$(EXEEXT)
- $(CXXLINK) $(pkveg2shadow_OBJECTS) $(pkveg2shadow_LDADD) $(LIBS)
mostlyclean-compile:
-rm -f *.$(OBJEXT)
@@ -534,6 +536,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pkcreatect.Po at am__quote@
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pkcrop.Po at am__quote@
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pkdiff.Po at am__quote@
+ at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pkdsm2shadow.Po at am__quote@
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pkdumpimg.Po at am__quote@
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pkdumpogr.Po at am__quote@
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pkegcs.Po at am__quote@
@@ -551,7 +554,6 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pksieve.Po at am__quote@
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pkstat.Po at am__quote@
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pkstatogr.Po at am__quote@
- at AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/pkveg2shadow.Po at am__quote@
@AMDEP_TRUE@@am__include@ @am__quote at ./$(DEPDIR)/svm.Po at am__quote@
.cc.o:
diff --git a/src/apps/pkclassify_svm.cc b/src/apps/pkclassify_svm.cc
new file mode 100644
index 0000000..25eb555
--- /dev/null
+++ b/src/apps/pkclassify_svm.cc
@@ -0,0 +1,1037 @@
+/**********************************************************************
+pkclassify_svm.cc: classify raster image using Artificial Neural Network
+Copyright (C) 2008-2012 Pieter Kempeneers
+
+This file is part of pktools
+
+pktools is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+pktools is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with pktools. If not, see <http://www.gnu.org/licenses/>.
+***********************************************************************/
+#include <vector>
+#include <map>
+#include <algorithm>
+#include "imageclasses/ImgReaderGdal.h"
+#include "imageclasses/ImgWriterGdal.h"
+#include "imageclasses/ImgReaderOgr.h"
+#include "imageclasses/ImgWriterOgr.h"
+#include "base/Optionpk.h"
+#include "algorithms/ConfusionMatrix.h"
+#include "algorithms/svm.h"
+#include "pkclassify_nn.h"
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#define Malloc(type,n) (type *)malloc((n)*sizeof(type))
+
+int main(int argc, char *argv[])
+{
+ map<short,int> reclassMap;
+ vector<int> vreclass;
+ vector<double> priors;
+
+ //--------------------------- command line options ------------------------------------
+
+ std::string versionString="version ";
+ versionString+=VERSION;
+ versionString+=", Copyright (C) 2008-2012 Pieter Kempeneers.\n\
+ This program comes with ABSOLUTELY NO WARRANTY; for details type use option -h.\n\
+ This is free software, and you are welcome to redistribute it\n\
+ under certain conditions; use option --license for details.";
+ Optionpk<bool> version_opt("\0","version",versionString,false);
+ Optionpk<bool> license_opt("lic","license","show license information",false);
+ Optionpk<bool> help_opt("h","help","shows this help info",false);
+ Optionpk<bool> todo_opt("\0","todo","",false);
+ Optionpk<string> input_opt("i", "input", "input image");
+ Optionpk<string> training_opt("t", "training", "training shape file. A single shape file contains all training features (must be set as: B0, B1, B2,...) for all classes (class numbers identified by label option). Use multiple training files for bootstrap aggregation (alternative to the bag and bsize options, where a random subset is taken from a single training file)");
+ Optionpk<string> label_opt("\0", "label", "identifier for class label in training shape file.","label");
+ Optionpk<unsigned short> reclass_opt("\0", "rc", "reclass code (e.g. --rc=12 --rc=23 to reclass first two classes to 12 and 23 resp.).", 0);
+ Optionpk<unsigned int> balance_opt("\0", "balance", "balance the input data to this number of samples for each class", 0);
+ Optionpk<int> minSize_opt("m", "min", "if number of training pixels is less then min, do not take this class into account", 0);
+ Optionpk<double> start_opt("s", "start", "start band sequence number (set to 0)",0);
+ Optionpk<double> end_opt("e", "end", "end band sequence number (set to 0 for all bands)", 0);
+ Optionpk<double> offset_opt("\0", "offset", "offset value for each spectral band input features: refl[band]=(DN[band]-offset[band])/scale[band]", 0.0);
+ Optionpk<double> scale_opt("\0", "scale", "scale value for each spectral band input features: refl=(DN[band]-offset[band])/scale[band] (use 0 if scale min and max in each band to -1.0 and 1.0)", 0.0);
+ Optionpk<unsigned short> aggreg_opt("a", "aggreg", "how to combine aggregated classifiers, see also rc option (0: no aggregation, 1: sum rule, 2: max rule).",0);
+ Optionpk<double> priors_opt("p", "prior", "prior probabilities for each class (e.g., -p 0.3 -p 0.3 -p 0.2 )", 0.0);
+
+
+ Optionpk<unsigned short> svm_type_opt("svmt", "svmtype", "type of SVM (0: C-SVC, 1: nu-SVC, 2: one-class SVM, 3: epsilon-SVR, 4: nu-SVR)",0);
+ Optionpk<unsigned short> kernel_type_opt("kt", "kerneltype", "type of kernel function (0: linear: u'*v, 1: polynomial: (gamma*u'*v + coef0)^degree, 2: radial basis function: exp(-gamma*|u-v|^2), 3: sigmoid: tanh(gamma*u'*v + coef0), 4: precomputed kernel (kernel values in training_set_file)",2);
+ Optionpk<unsigned short> kernel_degree_opt("kd", "kd", "degree in kernel function",3);
+ Optionpk<float> gamma_opt("g", "gamma", "gamma in kernel function",0);
+ Optionpk<float> coef0_opt("c0", "coef0", "coef0 in kernel function",0);
+ Optionpk<float> ccost_opt("cc", "ccost", "the parameter C of C-SVC, epsilon-SVR, and nu-SVR",1);
+ Optionpk<float> nu_opt("nu", "nu", "the parameter nu of nu-SVC, one-class SVM, and nu-SVR",0.5);
+ Optionpk<float> epsilon_loss_opt("eloss", "eloss", "the epsilon in loss function of epsilon-SVR",0.1);
+ Optionpk<int> cache_opt("cache", "cache", "cache memory size in MB",100);
+ Optionpk<float> epsilon_tol_opt("etol", "etol", "the tolerance of termination criterion",0.001);
+ Optionpk<bool> shrinking_opt("shrink", "shrink", "whether to use the shrinking heuristics",false);
+ Optionpk<bool> prob_est_opt("pe", "probest", "whether to train a SVC or SVR model for probability estimates",false);
+ // Optionpk<bool> weight_opt("wi", "wi", "set the parameter C of class i to weight*C, for C-SVC",true);
+ Optionpk<unsigned int> cv_opt("cv", "cv", "n-fold cross validation mode",0);
+ Optionpk<unsigned short> comb_opt("c", "comb", "how to combine bootstrap aggregation classifiers (0: sum rule, 1: product rule, 2: max rule). Also used to aggregate classes with rc option.",0);
+ Optionpk<unsigned short> bag_opt("\0", "bag", "Number of bootstrap aggregations", 1);
+ Optionpk<int> bagSize_opt("\0", "bsize", "Percentage of features used from available training features for each bootstrap aggregation", 100);
+ Optionpk<string> classBag_opt("\0", "class", "output for each individual bootstrap aggregation");
+ Optionpk<string> mask_opt("\0", "mask", "mask image (see also mvalue option");
+ Optionpk<short> maskValue_opt("\0", "mvalue", "mask value(s) not to consider for classification (use negative values if only these values should be taken into account). Values will be taken over in classification image.", 0);
+ Optionpk<unsigned short> flag_opt("f", "flag", "flag to put where image is invalid.", 0);
+ Optionpk<string> output_opt("o", "output", "output classification image");
+ Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate). Empty string: inherit from input image");
+ Optionpk<string> option_opt("co", "co", "options: NAME=VALUE [-co COMPRESS=LZW] [-co INTERLEAVE=BAND]");
+ Optionpk<string> colorTable_opt("\0", "ct", "colour table in ascii format having 5 columns: id R G B ALFA (0: transparent, 255: solid)");
+ Optionpk<string> prob_opt("\0", "prob", "probability image.");
+ Optionpk<short> verbose_opt("v", "verbose", "set to: 0 (results only), 1 (confusion matrix), 2 (debug)",0);
+
+ version_opt.retrieveOption(argc,argv);
+ license_opt.retrieveOption(argc,argv);
+ help_opt.retrieveOption(argc,argv);
+ todo_opt.retrieveOption(argc,argv);
+
+ input_opt.retrieveOption(argc,argv);
+ training_opt.retrieveOption(argc,argv);
+ label_opt.retrieveOption(argc,argv);
+ reclass_opt.retrieveOption(argc,argv);
+ balance_opt.retrieveOption(argc,argv);
+ minSize_opt.retrieveOption(argc,argv);
+ start_opt.retrieveOption(argc,argv);
+ end_opt.retrieveOption(argc,argv);
+ offset_opt.retrieveOption(argc,argv);
+ scale_opt.retrieveOption(argc,argv);
+ aggreg_opt.retrieveOption(argc,argv);
+ priors_opt.retrieveOption(argc,argv);
+ svm_type_opt.retrieveOption(argc,argv);
+ kernel_type_opt.retrieveOption(argc,argv);
+ kernel_degree_opt.retrieveOption(argc,argv);
+ gamma_opt.retrieveOption(argc,argv);
+ coef0_opt.retrieveOption(argc,argv);
+ ccost_opt.retrieveOption(argc,argv);
+ nu_opt.retrieveOption(argc,argv);
+ epsilon_loss_opt.retrieveOption(argc,argv);
+ cache_opt.retrieveOption(argc,argv);
+ epsilon_tol_opt.retrieveOption(argc,argv);
+ shrinking_opt.retrieveOption(argc,argv);
+ prob_est_opt.retrieveOption(argc,argv);
+ cv_opt.retrieveOption(argc,argv);
+ comb_opt.retrieveOption(argc,argv);
+ bag_opt.retrieveOption(argc,argv);
+ bagSize_opt.retrieveOption(argc,argv);
+ classBag_opt.retrieveOption(argc,argv);
+ mask_opt.retrieveOption(argc,argv);
+ maskValue_opt.retrieveOption(argc,argv);
+ flag_opt.retrieveOption(argc,argv);
+ output_opt.retrieveOption(argc,argv);
+ oformat_opt.retrieveOption(argc,argv);
+ colorTable_opt.retrieveOption(argc,argv);
+ option_opt.retrieveOption(argc,argv);
+ prob_opt.retrieveOption(argc,argv);
+ verbose_opt.retrieveOption(argc,argv);
+
+ if(version_opt[0]||todo_opt[0]){
+ std::cout << version_opt.getHelp() << std::endl;
+ std::cout << "todo: " << todo_opt.getHelp() << std::endl;
+ exit(0);
+ }
+ if(license_opt[0]){
+ std::cout << Optionpk<bool>::getGPLv3License() << std::endl;
+ exit(0);
+ }
+ if(help_opt[0]){
+ std::cout << "usage: pkclassify_nn -i testimage -o outputimage -t training [OPTIONS]" << std::endl;
+ exit(0);
+ }
+
+ if(verbose_opt[0]>=1){
+ std::cout << "image filename: " << input_opt[0] << std::endl;
+ if(mask_opt.size())
+ std::cout << "mask filename: " << mask_opt[0] << std::endl;
+ if(training_opt[0].size()){
+ std::cout << "training shape file: " << std::endl;
+ for(int ifile=0;ifile<training_opt.size();++ifile)
+ std::cout << training_opt[ifile] << std::endl;
+ }
+ else
+ cerr << "no training file set!" << std::endl;
+ std::cout << "verbose: " << verbose_opt[0] << std::endl;
+ }
+ unsigned short nbag=(training_opt.size()>1)?training_opt.size():bag_opt[0];
+ if(verbose_opt[0]>=1)
+ std::cout << "number of bootstrap aggregations: " << nbag << std::endl;
+
+ unsigned int totalSamples=0;
+ int nreclass=0;
+ vector<int> vcode;//unique class codes in recode string
+ // vector<FANN::neural_net> net(nbag);//the neural network
+ vector<struct svm_model*> svm(nbag);
+ vector<struct svm_parameter> param(nbag);
+
+ unsigned int nclass=0;
+ int nband=0;
+ int startBand=2;//first two bands represent X and Y pos
+
+ vector< vector<double> > offset(nbag);
+ vector< vector<double> > scale(nbag);
+ vector< Vector2d<float> > trainingPixels;//[class][sample][band]
+
+ if(reclass_opt.size()>1){
+ vreclass.resize(reclass_opt.size());
+ for(int iclass=0;iclass<reclass_opt.size();++iclass){
+ reclassMap[iclass]=reclass_opt[iclass];
+ vreclass[iclass]=reclass_opt[iclass];
+ }
+ }
+ if(priors_opt.size()>1){//priors from argument list
+ priors.resize(priors_opt.size());
+ double normPrior=0;
+ for(int iclass=0;iclass<priors_opt.size();++iclass){
+ priors[iclass]=priors_opt[iclass];
+ normPrior+=priors[iclass];
+ }
+ //normalize
+ for(int iclass=0;iclass<priors_opt.size();++iclass)
+ priors[iclass]/=normPrior;
+ }
+
+ //----------------------------------- Training -------------------------------
+ vector<struct svm_problem> prob(nbag);
+ vector<struct svm_node *> x_space(nbag);
+ // struct svm_node *x_space;
+ vector<string> fields;
+ for(int ibag=0;ibag<nbag;++ibag){
+ //organize training data
+ if(ibag<training_opt.size()){//if bag contains new training pixels
+ trainingPixels.clear();
+ map<int,Vector2d<float> > trainingMap;
+ if(verbose_opt[0]>=1)
+ std::cout << "reading imageShape file " << training_opt[0] << std::endl;
+ try{
+ totalSamples=readDataImageShape(training_opt[ibag],trainingMap,fields,start_opt[0],end_opt[0],label_opt[0],verbose_opt[0]);
+ if(trainingMap.size()<2){
+ string errorstring="Error: could not read at least two classes from training file";
+ throw(errorstring);
+ }
+ }
+ catch(string error){
+ cerr << error << std::endl;
+ exit(1);
+ }
+ catch(...){
+ cerr << "error catched" << std::endl;
+ exit(1);
+ }
+ //delete class 0
+ if(verbose_opt[0]>=1)
+ std::cout << "erasing class 0 from training set (" << trainingMap[0].size() << " from " << totalSamples << ") samples" << std::endl;
+ totalSamples-=trainingMap[0].size();
+ trainingMap.erase(0);
+ //convert map to vector
+ short iclass=0;
+ if(reclass_opt.size()==1){//no reclass option, read classes from shape
+ reclassMap.clear();
+ vreclass.clear();
+ }
+ if(verbose_opt[0]>1)
+ std::cout << "training pixels: " << std::endl;
+ map<int,Vector2d<float> >::iterator mapit=trainingMap.begin();
+ while(mapit!=trainingMap.end()){
+// for(map<int,Vector2d<float> >::const_iterator mapit=trainingMap.begin();mapit!=trainingMap.end();++mapit){
+ //delete small classes
+ if((mapit->second).size()<minSize_opt[0]){
+ trainingMap.erase(mapit);
+ continue;
+ //todo: beware of reclass option: delete this reclass if no samples are left in this classes!!
+ }
+ if(reclass_opt.size()==1){//no reclass option, read classes from shape
+ reclassMap[iclass]=(mapit->first);
+ vreclass.push_back(mapit->first);
+ }
+ trainingPixels.push_back(mapit->second);
+ if(verbose_opt[0]>1)
+ std::cout << mapit->first << ": " << (mapit->second).size() << " samples" << std::endl;
+ ++iclass;
+ ++mapit;
+ }
+ if(!ibag){
+ nclass=trainingPixels.size();
+ nband=trainingPixels[0][0].size()-2;//X and Y//trainingPixels[0][0].size();
+ }
+ else{
+ assert(nclass==trainingPixels.size());
+ assert(nband==trainingPixels[0][0].size()-2);
+ }
+ assert(reclassMap.size()==nclass);
+
+ //do not remove outliers here: could easily be obtained through ogr2ogr -where 'B2<110' output.shp input.shp
+ //balance training data
+ if(balance_opt[0]>0){
+ if(random)
+ srand(time(NULL));
+ totalSamples=0;
+ for(int iclass=0;iclass<nclass;++iclass){
+ if(trainingPixels[iclass].size()>balance_opt[0]){
+ while(trainingPixels[iclass].size()>balance_opt[0]){
+ int index=rand()%trainingPixels[iclass].size();
+ trainingPixels[iclass].erase(trainingPixels[iclass].begin()+index);
+ }
+ }
+ else{
+ int oldsize=trainingPixels[iclass].size();
+ for(int isample=trainingPixels[iclass].size();isample<balance_opt[0];++isample){
+ int index = rand()%oldsize;
+ trainingPixels[iclass].push_back(trainingPixels[iclass][index]);
+ }
+ }
+ totalSamples+=trainingPixels[iclass].size();
+ }
+ assert(totalSamples==nclass*balance_opt[0]);
+ }
+
+ //set scale and offset
+ offset[ibag].resize(nband);
+ scale[ibag].resize(nband);
+ if(offset_opt.size()>1)
+ assert(offset_opt.size()==nband);
+ if(scale_opt.size()>1)
+ assert(scale_opt.size()==nband);
+ Histogram hist;
+ for(int iband=0;iband<nband;++iband){
+ if(verbose_opt[0]>=1)
+ std::cout << "scaling for band" << iband << std::endl;
+ offset[ibag][iband]=(offset_opt.size()==1)?offset_opt[0]:offset_opt[iband];
+ scale[ibag][iband]=(scale_opt.size()==1)?scale_opt[0]:scale_opt[iband];
+ //search for min and maximum
+ if(scale[ibag][iband]<=0){
+ float theMin=trainingPixels[0][0][iband+startBand];
+ float theMax=trainingPixels[0][0][iband+startBand];
+ for(int iclass=0;iclass<nclass;++iclass){
+ for(int isample=0;isample<trainingPixels[iclass].size();++isample){
+ if(theMin>trainingPixels[iclass][isample][iband+startBand])
+ theMin=trainingPixels[iclass][isample][iband+startBand];
+ if(theMax<trainingPixels[iclass][isample][iband+startBand])
+ theMax=trainingPixels[iclass][isample][iband+startBand];
+ }
+ }
+ offset[ibag][iband]=theMin+(theMax-theMin)/2.0;
+ scale[ibag][iband]=(theMax-theMin)/2.0;
+ if(verbose_opt[0]>=1){
+ std::cout << "Extreme image values for band " << iband << ": [" << theMin << "," << theMax << "]" << std::endl;
+ std::cout << "Using offset, scale: " << offset[ibag][iband] << ", " << scale[ibag][iband] << std::endl;
+ std::cout << "scaled values for band " << iband << ": [" << (theMin-offset[ibag][iband])/scale[ibag][iband] << "," << (theMax-offset[ibag][iband])/scale[ibag][iband] << "]" << std::endl;
+ }
+ }
+ }
+ }
+ else{//use same offset and scale
+ offset[ibag].resize(nband);
+ scale[ibag].resize(nband);
+ for(int iband=0;iband<nband;++iband){
+ offset[ibag][iband]=offset[0][iband];
+ scale[ibag][iband]=scale[0][iband];
+ }
+ }
+
+ if(!ibag){
+ //recode vreclass to ordered vector, starting from 0 to nreclass
+ vcode.clear();
+ if(verbose_opt[0]>=1){
+ std::cout << "before recoding: " << std::endl;
+ for(int iclass = 0; iclass < vreclass.size(); iclass++)
+ std::cout << " " << vreclass[iclass];
+ std::cout << std::endl;
+ }
+ vector<int> vord=vreclass;//ordered vector, starting from 0 to nreclass
+ int iclass=0;
+ map<short,int> mreclass;
+ for(int ic=0;ic<vreclass.size();++ic){
+ if(mreclass.find(vreclass[ic])==mreclass.end())
+ mreclass[vreclass[ic]]=iclass++;
+ }
+ for(int ic=0;ic<vreclass.size();++ic)
+ vord[ic]=mreclass[vreclass[ic]];
+ //construct uniqe class codes
+ while(!vreclass.empty()){
+ vcode.push_back(*(vreclass.begin()));
+ //delete all these entries from vreclass
+ vector<int>::iterator vit;
+ while((vit=find(vreclass.begin(),vreclass.end(),vcode.back()))!=vreclass.end())
+ vreclass.erase(vit);
+ }
+ if(verbose_opt[0]>=1){
+ std::cout << "recode values: " << std::endl;
+ for(int icode=0;icode<vcode.size();++icode)
+ std::cout << vcode[icode] << " ";
+ std::cout << std::endl;
+ }
+ vreclass=vord;
+ if(verbose_opt[0]>=1){
+ std::cout << "after recoding: " << std::endl;
+ for(int iclass = 0; iclass < vord.size(); iclass++)
+ std::cout << " " << vord[iclass];
+ std::cout << std::endl;
+ }
+
+ vector<int> vuniqueclass=vreclass;
+ //remove duplicate elements from vuniqueclass
+ sort( vuniqueclass.begin(), vuniqueclass.end() );
+ vuniqueclass.erase( unique( vuniqueclass.begin(), vuniqueclass.end() ), vuniqueclass.end() );
+ nreclass=vuniqueclass.size();
+ if(verbose_opt[0]>=1){
+ std::cout << "unique classes: " << std::endl;
+ for(int iclass = 0; iclass < vuniqueclass.size(); iclass++)
+ std::cout << " " << vuniqueclass[iclass];
+ std::cout << std::endl;
+ std::cout << "number of reclasses: " << nreclass << std::endl;
+ }
+
+ if(priors_opt.size()==1){//default: equal priors for each class
+ priors.resize(nclass);
+ for(int iclass=0;iclass<nclass;++iclass)
+ priors[iclass]=1.0/nclass;
+ }
+ assert(priors_opt.size()==1||priors_opt.size()==nclass);
+
+ if(verbose_opt[0]>=1){
+ std::cout << "number of bands: " << nband << std::endl;
+ std::cout << "number of classes: " << nclass << std::endl;
+ std::cout << "priors:";
+ for(int iclass=0;iclass<nclass;++iclass)
+ std::cout << " " << priors[iclass];
+ std::cout << std::endl;
+ }
+ }//if(!ibag)
+
+ //Calculate features of trainig set
+ vector< Vector2d<float> > trainingFeatures(nclass);
+ for(int iclass=0;iclass<nclass;++iclass){
+ int nctraining=0;
+ if(verbose_opt[0]>=1)
+ std::cout << "calculating features for class " << iclass << std::endl;
+ if(random)
+ srand(time(NULL));
+ nctraining=(bagSize_opt[0]<100)? trainingPixels[iclass].size()/100.0*bagSize_opt[0] : trainingPixels[iclass].size();//bagSize_opt[0] given in % of training size
+ if(nctraining<=0)
+ nctraining=1;
+ assert(nctraining<=trainingPixels[iclass].size());
+ int index=0;
+ if(bagSize_opt[0]<100)
+ random_shuffle(trainingPixels[iclass].begin(),trainingPixels[iclass].end());
+
+ trainingFeatures[iclass].resize(nctraining);
+ for(int isample=0;isample<nctraining;++isample){
+ //scale pixel values according to scale and offset!!!
+ for(int iband=0;iband<nband;++iband){
+ float value=trainingPixels[iclass][isample][iband+startBand];
+ trainingFeatures[iclass][isample].push_back((value-offset[ibag][iband])/scale[ibag][iband]);
+ }
+ }
+ assert(trainingFeatures[iclass].size()==nctraining);
+ }
+
+ unsigned int nFeatures=trainingFeatures[0][0].size();
+ unsigned int ntraining=0;
+ for(int iclass=0;iclass<nclass;++iclass){
+ if(verbose_opt[0]>=1)
+ std::cout << "training sample size for class " << vcode[iclass] << ": " << trainingFeatures[iclass].size() << std::endl;
+ ntraining+=trainingFeatures[iclass].size();
+ }
+ // vector<struct svm_problem> prob(ibag);
+ // vector<struct svm_node *> x_space(ibag);
+ prob[ibag].l=ntraining;
+ prob[ibag].y = Malloc(double,prob[ibag].l);
+ prob[ibag].x = Malloc(struct svm_node *,prob[ibag].l);
+ x_space[ibag] = Malloc(struct svm_node,(nFeatures+1)*ntraining);
+ unsigned long int spaceIndex=0;
+ int lIndex=0;
+ for(int iclass=0;iclass<nclass;++iclass){
+ for(int isample=0;isample<trainingFeatures[iclass].size();++isample){
+ // //test
+ // std::cout << iclass;
+ prob[ibag].x[lIndex]=&(x_space[ibag][spaceIndex]);
+ for(int ifeature=0;ifeature<nFeatures;++ifeature){
+ x_space[ibag][spaceIndex].index=ifeature+1;
+ x_space[ibag][spaceIndex].value=trainingFeatures[iclass][isample][ifeature];
+ // //test
+ // std::cout << " " << x_space[ibag][spaceIndex].index << ":" << x_space[ibag][spaceIndex].value;
+ ++spaceIndex;
+ }
+ x_space[ibag][spaceIndex++].index=-1;
+ prob[ibag].y[lIndex]=iclass;
+ ++lIndex;
+ }
+ }
+ assert(lIndex==prob[ibag].l);
+
+ //set SVM parameters through command line options
+ param[ibag].svm_type = svm_type_opt[0];
+ param[ibag].kernel_type = kernel_type_opt[0];
+ param[ibag].degree = kernel_degree_opt[0];
+ param[ibag].gamma = (gamma_opt[0]>0)? gamma_opt[0] : 1.0/nFeatures;
+ param[ibag].coef0 = coef0_opt[0];
+ param[ibag].nu = nu_opt[0];
+ param[ibag].cache_size = cache_opt[0];
+ param[ibag].C = ccost_opt[0];
+ param[ibag].eps = epsilon_tol_opt[0];
+ param[ibag].p = epsilon_loss_opt[0];
+ param[ibag].shrinking = (shrinking_opt[0])? 1 : 0;
+ param[ibag].probability = (prob_est_opt[0])? 1 : 0;
+ param[ibag].nr_weight = 0;//not used: I use priors and balancing
+ param[ibag].weight_label = NULL;
+ param[ibag].weight = NULL;
+
+ if(verbose_opt[0])
+ std::cout << "checking parameters" << std::endl;
+ svm_check_parameter(&prob[ibag],¶m[ibag]);
+ if(verbose_opt[0])
+ std::cout << "parameters ok, training" << std::endl;
+ svm[ibag]=svm_train(&prob[ibag],¶m[ibag]);
+
+ if(verbose_opt[0])
+ std::cout << "SVM is now trained" << std::endl;
+ if(cv_opt[0]>0){
+ std::cout << "Confusion matrix" << std::endl;
+ ConfusionMatrix cm(nclass);
+ // for(int iclass=0;iclass<nclass;++iclass)
+ // cm.pushBackClassName(type2string(iclass));
+ double *target = Malloc(double,prob[ibag].l);
+ svm_cross_validation(&prob[ibag],¶m[ibag],cv_opt[0],target);
+ assert(param[ibag].svm_type != EPSILON_SVR&¶m[ibag].svm_type != NU_SVR);//only for regression
+ int total_correct=0;
+ for(int i=0;i<prob[ibag].l;i++)
+ cm.incrementResult(cm.getClass(prob[ibag].y[i]),cm.getClass(target[i]),1);
+ assert(cm.nReference());
+ std::cout << cm << std::endl;
+ std::cout << "Kappa: " << cm.kappa() << std::endl;
+ double se95_oa=0;
+ double doa=0;
+ doa=cm.oa_pct(&se95_oa);
+ std::cout << "Overall Accuracy: " << doa << " (" << se95_oa << ")" << std::endl;
+ free(target);
+ }
+
+ // *NOTE* Because svm_model contains pointers to svm_problem, you can
+ // not free the memory used by svm_problem if you are still using the
+ // svm_model produced by svm_train().
+
+ // free(prob.y);
+ // free(prob.x);
+ // free(x_space);
+ // svm_destroy_param(¶m);
+ }//for ibag
+
+ //--------------------------------- end of training -----------------------------------
+ if(!output_opt.size())
+ exit(0);
+
+
+ const char* pszMessage;
+ void* pProgressArg=NULL;
+ GDALProgressFunc pfnProgress=GDALTermProgress;
+ float progress=0;
+ if(!verbose_opt[0])
+ pfnProgress(progress,pszMessage,pProgressArg);
+ //-------------------------------- open image file ------------------------------------
+ if(input_opt[0].find(".shp")==string::npos){
+ ImgReaderGdal testImage;
+ try{
+ if(verbose_opt[0]>=1)
+ std::cout << "opening image " << input_opt[0] << std::endl;
+ testImage.open(input_opt[0]);
+ }
+ catch(string error){
+ cerr << error << std::endl;
+ exit(2);
+ }
+ ImgReaderGdal maskReader;
+ if(mask_opt.size()){
+ try{
+ if(verbose_opt[0]>=1)
+ std::cout << "opening mask image file " << mask_opt[0] << std::endl;
+ maskReader.open(mask_opt[0]);
+ assert(maskReader.nrOfCol()==testImage.nrOfCol());
+ assert(maskReader.nrOfRow()==testImage.nrOfRow());
+ }
+ catch(string error){
+ cerr << error << std::endl;
+ exit(2);
+ }
+ catch(...){
+ cerr << "error catched" << std::endl;
+ exit(1);
+ }
+ }
+ int nrow=testImage.nrOfRow();
+ int ncol=testImage.nrOfCol();
+ if(option_opt.findSubstring("INTERLEAVE=")==option_opt.end()){
+ string theInterleave="INTERLEAVE=";
+ theInterleave+=testImage.getInterleave();
+ option_opt.push_back(theInterleave);
+ }
+ vector<char> classOut(ncol);//classified line for writing to image file
+
+ // assert(nband==testImage.nrOfBand());
+ ImgWriterGdal classImageBag;
+ ImgWriterGdal classImageOut;
+ ImgWriterGdal probImage;
+ string imageType=testImage.getImageType();
+ if(oformat_opt.size())//default
+ imageType=oformat_opt[0];
+ try{
+ assert(output_opt.size());
+ if(verbose_opt[0]>=1)
+ std::cout << "opening class image for writing output " << output_opt[0] << std::endl;
+ if(classBag_opt.size()){
+ classImageBag.open(output_opt[0],ncol,nrow,nbag,GDT_Byte,imageType,option_opt);
+ classImageBag.copyGeoTransform(testImage);
+ classImageBag.setProjection(testImage.getProjection());
+ }
+ classImageOut.open(output_opt[0],ncol,nrow,1,GDT_Byte,imageType,option_opt);
+ classImageOut.copyGeoTransform(testImage);
+ classImageOut.setProjection(testImage.getProjection());
+ if(colorTable_opt.size())
+ classImageOut.setColorTable(colorTable_opt[0],0);
+ if(prob_opt.size()){
+ probImage.open(prob_opt[0],ncol,nrow,nreclass,GDT_Byte,imageType,option_opt);
+ probImage.copyGeoTransform(testImage);
+ probImage.setProjection(testImage.getProjection());
+ }
+ }
+ catch(string error){
+ cerr << error << std::endl;
+ }
+
+ for(int iline=0;iline<nrow;++iline){
+ vector<float> buffer(ncol);
+ vector<short> lineMask;
+ if(mask_opt.size())
+ lineMask.resize(maskReader.nrOfCol());
+ Vector2d<float> hpixel(ncol,nband);
+ // Vector2d<float> fpixel(ncol);
+ Vector2d<float> prOut(nreclass,ncol);//posterior prob for each reclass
+ Vector2d<char> classBag;//classified line for writing to image file
+ if(classBag_opt.size())
+ classBag.resize(nbag,ncol);
+ //read all bands of all pixels in this line in hline
+ for(int iband=start_opt[0];iband<start_opt[0]+nband;++iband){
+ if(verbose_opt[0]==2)
+ std::cout << "reading band " << iband << std::endl;
+ assert(iband>=0);
+ assert(iband<testImage.nrOfBand());
+ try{
+ testImage.readData(buffer,GDT_Float32,iline,iband);
+ }
+ catch(string theError){
+ cerr << "Error reading " << input_opt[0] << ": " << theError << std::endl;
+ exit(3);
+ }
+ catch(...){
+ cerr << "error catched" << std::endl;
+ exit(3);
+ }
+ for(int icol=0;icol<ncol;++icol)
+ hpixel[icol][iband-start_opt[0]]=buffer[icol];
+ }
+
+ assert(nband==hpixel[0].size());
+ if(verbose_opt[0]==2)
+ std::cout << "used bands: " << nband << std::endl;
+ //read mask
+ if(!lineMask.empty()){
+ try{
+ maskReader.readData(lineMask,GDT_Int16,iline);
+ }
+ catch(string theError){
+ cerr << "Error reading " << mask_opt[0] << ": " << theError << std::endl;
+ exit(3);
+ }
+ catch(...){
+ cerr << "error catched" << std::endl;
+ exit(3);
+ }
+ }
+
+ //process per pixel
+ for(int icol=0;icol<ncol;++icol){
+
+ bool masked=false;
+ if(!lineMask.empty()){
+ short theMask=0;
+ for(short ivalue=0;ivalue<maskValue_opt.size();++ivalue){
+ if(maskValue_opt[ivalue]>=0){//values set in maskValue_opt are invalid
+ if(lineMask[icol]==maskValue_opt[ivalue]){
+ theMask=(flag_opt.size()==maskValue_opt.size())? flag_opt[ivalue] : flag_opt[0];// lineMask[icol];
+ masked=true;
+ break;
+ }
+ }
+ else{//only values set in maskValue_opt are valid
+ if(lineMask[icol]!=-maskValue_opt[ivalue]){
+ theMask=(flag_opt.size()==maskValue_opt.size())? flag_opt[ivalue] : flag_opt[0];// lineMask[icol];
+ masked=true;
+ }
+ else{
+ masked=false;
+ break;
+ }
+ }
+ }
+ if(masked){
+ if(classBag_opt.size())
+ for(int ibag=0;ibag<nbag;++ibag)
+ classBag[ibag][icol]=theMask;
+ classOut[icol]=theMask;
+ continue;
+ }
+ }
+ bool valid=false;
+ for(int iband=0;iband<nband;++iband){
+ if(hpixel[icol][iband]){
+ valid=true;
+ break;
+ }
+ }
+ if(!valid){
+ if(classBag_opt.size())
+ for(int ibag=0;ibag<nbag;++ibag)
+ classBag[ibag][icol]=flag_opt[0];
+ classOut[icol]=flag_opt[0];
+ continue;//next column
+ }
+ for(int iclass=0;iclass<nreclass;++iclass)
+ prOut[iclass][icol]=0;
+ //----------------------------------- classification -------------------
+ for(int ibag=0;ibag<nbag;++ibag){
+ //calculate image features
+ // fpixel[icol].clear();
+ // for(int iband=0;iband<nband;++iband)
+ // fpixel[icol].push_back((hpixel[icol][iband]-offset[ibag][iband])/scale[ibag][iband]);
+ vector<double> result(nclass);
+ // result=net[ibag].run(fpixel[icol]);
+ struct svm_node *x;
+ // x = (struct svm_node *) malloc((fpixel[icol].size()+1)*sizeof(struct svm_node));
+ x = (struct svm_node *) malloc((nband+1)*sizeof(struct svm_node));
+ // for(int i=0;i<fpixel[icol].size();++i){
+ for(int iband=0;iband<nband;++iband){
+ x[iband].index=iband+1;
+ // x[i].value=fpixel[icol][i];
+ x[iband].value=(hpixel[icol][iband]-offset[ibag][iband])/scale[ibag][iband];
+ }
+ // x[fpixel[icol].size()].index=-1;//to end svm feature vector
+ x[nband].index=-1;//to end svm feature vector
+ double predict_label=0;
+ vector<float> pValues(nclass);
+ vector<float> prValues(nreclass);
+ vector<float> priorsReclass(nreclass);
+ float maxP=0;
+ if(!aggreg_opt[0]){
+ predict_label = svm_predict(svm[ibag],x);
+ for(int iclass=0;iclass<nclass;++iclass){
+ if(iclass==static_cast<int>(predict_label))
+ result[iclass]=1;
+ else
+ result[iclass]=0;
+ }
+ }
+ else{
+ assert(svm_check_probability_model(svm[ibag]));
+ predict_label = svm_predict_probability(svm[ibag],x,&(result[0]));
+ }
+ for(int iclass=0;iclass<nclass;++iclass){
+ float pv=result[iclass];
+ assert(pv>=0);
+ assert(pv<=1);
+ pv*=priors[iclass];
+ pValues[iclass]=pv;
+ }
+ float normReclass=0;
+ for(int iclass=0;iclass<nreclass;++iclass){
+ prValues[iclass]=0;
+ priorsReclass[iclass]=0;
+ float maxPaggreg=0;
+ for(int ic=0;ic<nclass;++ic){
+ if(vreclass[ic]==iclass){
+ priorsReclass[iclass]+=priors[ic];
+ switch(aggreg_opt[0]){
+ default:
+ case(1)://sum rule (sum posterior probabilities of aggregated individual classes)
+ prValues[iclass]+=pValues[ic];
+ break;
+ case(0):
+ case(2)://max rule (look for maximum post probability of aggregated individual classes)
+ if(pValues[ic]>maxPaggreg){
+ maxPaggreg=pValues[ic];
+ prValues[iclass]=maxPaggreg;
+ }
+ break;
+ }
+ }
+ }
+ }
+ for(int iclass=0;iclass<nreclass;++iclass)
+ normReclass+=prValues[iclass];
+
+ //calculate posterior prob of bag
+ if(classBag_opt.size()){
+ //search for max prob within bag
+ maxP=0;
+ classBag[ibag][icol]=0;
+ }
+ for(int iclass=0;iclass<nreclass;++iclass){
+ float prv=prValues[iclass];
+ prv/=normReclass;
+ // prv*=100.0;
+ prValues[iclass]=prv;
+ switch(comb_opt[0]){
+ default:
+ case(0)://sum rule
+ prOut[iclass][icol]+=prValues[iclass]+static_cast<float>(1.0-nbag)/nbag*priorsReclass[iclass];//add probabilities for each bag
+ break;
+ case(1)://product rule
+ prOut[iclass][icol]*=pow(priorsReclass[iclass],static_cast<float>(1.0-nbag)/nbag)*prValues[iclass];//add probabilities for each bag
+ break;
+ case(2)://max rule
+ if(prValues[iclass]>prOut[iclass][icol])
+ prOut[iclass][icol]=prValues[iclass];
+ break;
+ }
+ if(classBag_opt.size()){
+ //search for max prob within bag
+ if(prValues[iclass]>maxP){
+ maxP=prValues[iclass];
+ classBag[ibag][icol]=vcode[iclass];
+ }
+ }
+ }
+ free(x);
+ }//ibag
+
+ //search for max class prob
+ float maxBag=0;
+ float normBag=0;
+ for(int iclass=0;iclass<nreclass;++iclass){
+ if(prOut[iclass][icol]>maxBag){
+ maxBag=prOut[iclass][icol];
+ classOut[icol]=vcode[iclass];
+ }
+ normBag+=prOut[iclass][icol];
+ }
+ //normalize prOut and convert to percentage
+ if(prob_opt.size()){
+ for(int iclass=0;iclass<nreclass;++iclass){
+ float prv=prOut[iclass][icol];
+ prv/=normBag;
+ prv*=100.0;
+ prOut[iclass][icol]=static_cast<short>(prv+0.5);
+ }
+ }
+ }//icol
+ //----------------------------------- write output ------------------------------------------
+ if(classBag_opt.size())
+ for(int ibag=0;ibag<nbag;++ibag)
+ classImageBag.writeData(classBag[ibag],GDT_Byte,iline,ibag);
+ if(prob_opt.size()){
+ for(int iclass=0;iclass<nreclass;++iclass)
+ probImage.writeData(prOut[iclass],GDT_Float32,iline,iclass);
+ }
+ classImageOut.writeData(classOut,GDT_Byte,iline);
+ if(!verbose_opt[0]){
+ progress=static_cast<float>(iline+1.0)/classImageOut.nrOfRow();
+ pfnProgress(progress,pszMessage,pProgressArg);
+ }
+ }
+ testImage.close();
+ if(prob_opt.size())
+ probImage.close();
+ if(classBag_opt.size())
+ classImageBag.close();
+ classImageOut.close();
+ }
+ else{//classify shape file
+ for(int ivalidation=0;ivalidation<input_opt.size();++ivalidation){
+ assert(output_opt.size()==input_opt.size());
+ if(verbose_opt[0])
+ std::cout << "opening img reader " << input_opt[ivalidation] << std::endl;
+ ImgReaderOgr imgReaderOgr(input_opt[ivalidation]);
+ if(verbose_opt[0])
+ std::cout << "opening img writer and copying fields from img reader" << output_opt[ivalidation] << std::endl;
+ ImgWriterOgr imgWriterOgr(output_opt[ivalidation],imgReaderOgr,false);
+ if(verbose_opt[0])
+ std::cout << "creating field class" << std::endl;
+ imgWriterOgr.createField("class",OFTInteger);
+ OGRFeature *poFeature;
+ unsigned int ifeature=0;
+ unsigned int nFeatures=imgReaderOgr.getFeatureCount();
+ while( (poFeature = imgReaderOgr.getLayer()->GetNextFeature()) != NULL ){
+ if(verbose_opt[0]>1)
+ std::cout << "feature " << ifeature << std::endl;
+ if( poFeature == NULL )
+ break;
+ OGRFeature *poDstFeature = NULL;
+ poDstFeature=imgWriterOgr.createFeature();
+ if( poDstFeature->SetFrom( poFeature, TRUE ) != OGRERR_NONE ){
+ CPLError( CE_Failure, CPLE_AppDefined,
+ "Unable to translate feature %d from layer %s.\n",
+ poFeature->GetFID(), imgWriterOgr.getLayerName().c_str() );
+ OGRFeature::DestroyFeature( poFeature );
+ OGRFeature::DestroyFeature( poDstFeature );
+ }
+ vector<float> validationPixel;
+ vector<float> validationFeature;
+
+ imgReaderOgr.readData(validationPixel,OFTReal,fields,poFeature);
+ OGRFeature::DestroyFeature( poFeature );
+// assert(validationPixel.size()>=start_opt[0]+nband);
+ assert(validationPixel.size()==nband);
+ vector<float> prOut(nreclass);//posterior prob for each reclass
+ for(int iclass=0;iclass<nreclass;++iclass)
+ prOut[iclass]=0;
+ for(int ibag=0;ibag<nbag;++ibag){
+// for(int iband=start_opt[0];iband<start_opt[0]+nband;++iband){
+ for(int iband=0;iband<nband;++iband){
+// validationFeature.push_back((validationPixel[iband]-offset[ibag][iband-start_opt[0]])/scale[ibag][iband-start_opt[0]]);
+ validationFeature.push_back((validationPixel[iband]-offset[ibag][iband])/scale[ibag][iband]);
+ if(verbose_opt[0]==2)
+ std::cout << " " << validationFeature.back();
+ }
+ if(verbose_opt[0]==2)
+ std::cout << std::endl;
+ vector<double> result(nclass);
+
+ // result=net[ibag].run(validationFeature);
+ struct svm_node *x;
+ x = (struct svm_node *) malloc((validationFeature.size()+1)*sizeof(struct svm_node));
+ for(int i=0;i<validationFeature.size();++i){
+ x[i].index=i+1;
+ x[i].value=validationFeature[i];
+ }
+ x[validationFeature.size()].index=-1;//to end svm feature vector
+ double predict_label=0;
+ vector<float> pValues(nclass);
+ vector<float> prValues(nreclass);
+ vector<float> priorsReclass(nreclass);
+ float maxP=0;
+ if(!aggreg_opt[0]){
+ predict_label = svm_predict(svm[ibag],x);
+ for(int iclass=0;iclass<nclass;++iclass){
+ if(predict_label==static_cast<int>(predict_label))
+ result[iclass]=1;
+ else
+ result[iclass]=0;
+ }
+ }
+ else{
+ assert(svm_check_probability_model(svm[ibag]));
+ predict_label = svm_predict_probability(svm[ibag],x,&(result[0]));
+ }
+ // int maxClass=0;
+ for(int iclass=0;iclass<nclass;++iclass){
+ float pv=(result[iclass]+1.0)/2.0;//bring back to scale [0,1]
+ pv*=priors[iclass];
+ pValues[iclass]=pv;
+ }
+ float normReclass=0;
+ for(int iclass=0;iclass<nreclass;++iclass){
+ prValues[iclass]=0;
+ priorsReclass[iclass]=0;
+ float maxPaggreg=0;
+ for(int ic=0;ic<nclass;++ic){
+ if(vreclass[ic]==iclass){
+ priorsReclass[iclass]+=priors[ic];
+ switch(aggreg_opt[0]){
+ default:
+ case(0)://sum rule (sum posterior probabilities of aggregated individual classes)
+ prValues[iclass]+=pValues[ic];
+ break;
+ case(1)://max rule (look for maximum post probability of aggregated individual classes)
+ if(pValues[ic]>maxPaggreg){
+ maxPaggreg=pValues[ic];
+ prValues[iclass]=maxPaggreg;
+ }
+ break;
+ }
+ }
+ }
+ }
+ for(int iclass=0;iclass<nreclass;++iclass)
+ normReclass+=prValues[iclass];
+ //calculate posterior prob of bag
+ for(int iclass=0;iclass<nreclass;++iclass){
+ float prv=prValues[iclass];
+ prv/=normReclass;
+ // prv*=100.0;
+ prValues[iclass]=prv;
+ switch(comb_opt[0]){
+ default:
+ case(0)://sum rule
+ prOut[iclass]+=prValues[iclass]+static_cast<float>(1.0-nbag)/nbag*priorsReclass[iclass];//add probabilities for each bag
+ break;
+ case(1)://product rule
+ prOut[iclass]*=pow(priorsReclass[iclass],static_cast<float>(1.0-nbag)/nbag)*prValues[iclass];//add probabilities for each bag
+ break;
+ case(2)://max rule
+ if(prValues[iclass]>prOut[iclass])
+ prOut[iclass]=prValues[iclass];
+ break;
+ }
+ }
+ free(x);
+ }//for ibag
+ //search for max class prob
+ float maxBag=0;
+ float normBag=0;
+ char classOut=0;
+ for(int iclass=0;iclass<nreclass;++iclass){
+ if(prOut[iclass]>maxBag){
+ maxBag=prOut[iclass];
+ classOut=vcode[iclass];
+ }
+ normBag+=prOut[iclass];
+ }
+ //normalize prOut and convert to percentage
+ for(int iclass=0;iclass<nreclass;++iclass){
+ float prv=prOut[iclass];
+ prv/=normBag;
+ prv*=100.0;
+ prOut[iclass]=static_cast<short>(prv+0.5);
+ }
+ poDstFeature->SetField("class",classOut);
+ poDstFeature->SetFID( poFeature->GetFID() );
+ CPLErrorReset();
+ if(imgWriterOgr.createFeature( poDstFeature ) != OGRERR_NONE){
+ CPLError( CE_Failure, CPLE_AppDefined,
+ "Unable to translate feature %d from layer %s.\n",
+ poFeature->GetFID(), imgWriterOgr.getLayerName().c_str() );
+ OGRFeature::DestroyFeature( poDstFeature );
+ OGRFeature::DestroyFeature( poDstFeature );
+ }
+ ++ifeature;
+ if(!verbose_opt[0]){
+ progress=static_cast<float>(ifeature+1.0)/nFeatures;
+ pfnProgress(progress,pszMessage,pProgressArg);
+ }
+ }
+ imgReaderOgr.close();
+ imgWriterOgr.close();
+ }
+ }
+
+ for(int ibag=0;ibag<nbag;++ibag){
+ svm_destroy_param[ibag](¶m[ibag]);
+ free(prob[ibag].y);
+ free(prob[ibag].x);
+ free(x_space[ibag]);
+ svm_free_and_destroy_model(&(svm[ibag]));
+ }
+ return 0;
+}
diff --git a/src/apps/pkdsm2shadow.cc b/src/apps/pkdsm2shadow.cc
new file mode 100644
index 0000000..87dbd6f
--- /dev/null
+++ b/src/apps/pkdsm2shadow.cc
@@ -0,0 +1,148 @@
+/**********************************************************************
+pkveg2shadow.cc: program to calculate sun shadow based on digital surface model and sun angles)
+Copyright (C) 2008-2012 Pieter Kempeneers
+
+This file is part of pktools
+
+pktools is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+pktools is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with pktools. If not, see <http://www.gnu.org/licenses/>.
+***********************************************************************/
+#include <assert.h>
+#include <iostream>
+#include <string>
+#include <fstream>
+#include <math.h>
+#include <sys/types.h>
+#include <stdio.h>
+#include "base/Optionpk.h"
+#include "base/Vector2d.h"
+#include "algorithms/Filter2d.h"
+#include "imageclasses/ImgReaderGdal.h"
+#include "imageclasses/ImgWriterGdal.h"
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+/*------------------
+ Main procedure
+ ----------------*/
+int main(int argc,char **argv) {
+ std::string versionString="version ";
+ versionString+=VERSION;
+ versionString+=", Copyright (C) 2008-2012 Pieter Kempeneers.\n\
+ This program comes with ABSOLUTELY NO WARRANTY; for details type use option -h.\n\
+ This is free software, and you are welcome to redistribute it\n\
+ under certain conditions; use option --license for details.";
+ Optionpk<bool> version_opt("\0","version",versionString,false);
+ Optionpk<bool> license_opt("lic","license","show license information",false);
+ Optionpk<bool> help_opt("h","help","shows this help info",false);
+ Optionpk<bool> todo_opt("\0","todo","",false);
+ Optionpk<std::string> input_opt("i","input","input image file","");
+ Optionpk<std::string> output_opt("o", "output", "Output image file", "");
+ Optionpk<double> sza_opt("sza", "sza", "Sun zenith angle.");
+ Optionpk<double> saa_opt("saa", "saa", "Sun azimuth angle (N=0 E=90 S=180 W=270).");
+ Optionpk<int> flag_opt("f", "flag", "Flag to put in image if pixel shadow", 0);
+ Optionpk<std::string> otype_opt("ot", "otype", "Data type for output image ({Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/CInt16/CInt32/CFloat32/CFloat64}). Empty string: inherit type from input image", "");
+ Optionpk<string> oformat_opt("of", "oformat", "Output image format (see also gdal_translate). Empty string: inherit from input image");
+ Optionpk<string> colorTable_opt("ct", "ct", "color table (file with 5 columns: id R G B ALFA (0: transparent, 255: solid)", "");
+ Optionpk<std::string> option_opt("co", "co", "options: NAME=VALUE [-co COMPRESS=LZW] [-co INTERLEAVE=BAND]");
+ Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0);
+
+ version_opt.retrieveOption(argc,argv);
+ license_opt.retrieveOption(argc,argv);
+ help_opt.retrieveOption(argc,argv);
+ todo_opt.retrieveOption(argc,argv);
+
+ input_opt.retrieveOption(argc,argv);
+ output_opt.retrieveOption(argc,argv);
+ sza_opt.retrieveOption(argc,argv);
+ saa_opt.retrieveOption(argc,argv);
+ flag_opt.retrieveOption(argc,argv);
+ option_opt.retrieveOption(argc,argv);
+ otype_opt.retrieveOption(argc,argv);
+ oformat_opt.retrieveOption(argc,argv);
+ colorTable_opt.retrieveOption(argc,argv);
+ verbose_opt.retrieveOption(argc,argv);
+
+ if(version_opt[0]||todo_opt[0]){
+ cout << version_opt.getHelp() << endl;
+ cout << "todo: " << todo_opt.getHelp() << endl;
+ exit(0);
+ }
+ if(license_opt[0]){
+ cout << Optionpk<bool>::getGPLv3License() << endl;
+ exit(0);
+ }
+ if(help_opt[0]){
+ cout << "usage: pkveg2shadow -i inputimage -o outputimage [OPTIONS]" << endl;
+ exit(0);
+ }
+
+ ImgReaderGdal input;
+ ImgWriterGdal output;
+ input.open(input_opt[0]);
+ // output.open(output_opt[0],input);
+ GDALDataType theType=GDT_Unknown;
+ if(verbose_opt[0])
+ cout << "possible output data types: ";
+ for(int iType = 0; iType < GDT_TypeCount; ++iType){
+ if(verbose_opt[0])
+ cout << " " << GDALGetDataTypeName((GDALDataType)iType);
+ if( GDALGetDataTypeName((GDALDataType)iType) != NULL
+ && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
+ otype_opt[0].c_str()))
+ theType=(GDALDataType) iType;
+ }
+ if(theType==GDT_Unknown)
+ theType=input.getDataType();
+
+ if(verbose_opt[0])
+ std::cout << std::endl << "Output pixel type: " << GDALGetDataTypeName(theType) << endl;
+
+ string imageType=input.getImageType();
+ if(oformat_opt.size())
+ imageType=oformat_opt[0];
+
+ if(option_opt.findSubstring("INTERLEAVE=")==option_opt.end()){
+ string theInterleave="INTERLEAVE=";
+ theInterleave+=input.getInterleave();
+ option_opt.push_back(theInterleave);
+ }
+ try{
+ output.open(output_opt[0],input.nrOfCol(),input.nrOfRow(),input.nrOfBand(),theType,imageType,option_opt);
+ }
+ catch(string errorstring){
+ cout << errorstring << endl;
+ exit(4);
+ }
+ if(input.isGeoRef()){
+ output.setProjection(input.getProjection());
+ double ulx,uly,deltaX,deltaY,rot1,rot2;
+ input.getGeoTransform(ulx,uly,deltaX,deltaY,rot1,rot2);
+ output.setGeoTransform(ulx,uly,deltaX,deltaY,rot1,rot2);
+ }
+ if(input.getColorTable()!=NULL)
+ output.setColorTable(input.getColorTable());
+
+ Filter2d::Filter2d filter2d;
+ if(verbose_opt[0])
+ std::cout<< "class values: ";
+ if(colorTable_opt[0]!="")
+ output.setColorTable(colorTable_opt[0]);
+ filter2d.shadowDsm(input,output,sza_opt[0],saa_opt[0],input.getDeltaX(),flag_opt[0]);
+ input.close();
+ output.close();
+ return 0;
+}
diff --git a/src/apps/pkextract.cc b/src/apps/pkextract.cc
index b58ca44..8f63472 100644
--- a/src/apps/pkextract.cc
+++ b/src/apps/pkextract.cc
@@ -70,7 +70,7 @@ int main(int argc, char *argv[])
Optionpk<string> ltype_opt("lt", "ltype", "Label type: In16 or String", "Integer");
Optionpk<string> fieldname_opt("bn", "bname", "Field name of output shape file", "B");
Optionpk<string> label_opt("cn", "cname", "name of the class label in the output vector file", "label");
- Optionpk<bool> polygon_opt("l", "line", "create OGRPolygon as geometry instead of points. Only if sample option is also of polygon type. Use 0 for OGRPoint", 0);
+ Optionpk<bool> polygon_opt("l", "line", "create OGRPolygon as geometry instead of points. Only if sample option is also of polygon type. Use 0 for OGRPoint", false);
Optionpk<int> band_opt("b", "band", "band index to crop. Use -1 to use all bands)", -1);
Optionpk<short> rule_opt("r", "rule", "rule how to report image information per feature. 0: value at each point (or at centroid of the polygon if line is not set), 1: mean value (written to centroid of polygon if line is not set), 2: proportion classes, 3: custom, 4: minimum of polygon).", 0);
Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0);
@@ -481,7 +481,7 @@ int main(int argc, char *argv[])
std::cout << "class " << class_opt[iclass] << " has " << nvalid[iclass] << " samples" << std::endl;
}
}
- else{//classification file
+ else{//class_opt[0]!=0
assert(class_opt[0]);
// if(class_opt[0]){
assert(threshold_opt.size()==1||threshold_opt.size()==class_opt.size());
diff --git a/src/imageclasses/ImgWriterGdal.cc b/src/imageclasses/ImgWriterGdal.cc
index bbbdacf..06e9cb2 100644
--- a/src/imageclasses/ImgWriterGdal.cc
+++ b/src/imageclasses/ImgWriterGdal.cc
@@ -106,6 +106,7 @@ void ImgWriterGdal::setCodec(const ImgReaderGdal& imgSrc){
}
char **papszMetadata;
papszMetadata = poDriver->GetMetadata();
+ //todo: try and catch if CREATE is not supported (as in PNG)
assert( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATE, FALSE ));
char **papszOptions=NULL;
for(vector<string>::const_iterator optionIt=m_options.begin();optionIt!=m_options.end();++optionIt)
@@ -181,6 +182,7 @@ void ImgWriterGdal::setCodec(const string& imageType)
}
char **papszMetadata;
papszMetadata = poDriver->GetMetadata();
+ //todo: try and catch if CREATE is not supported (as in PNG)
assert( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATE, FALSE ));
char **papszOptions=NULL;
for(vector<string>::const_iterator optionIt=m_options.begin();optionIt!=m_options.end();++optionIt)
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-grass/pktools.git
More information about the Pkg-grass-devel
mailing list