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 <kempe...@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@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.@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@./$(DEPDIR)/pkcreatect.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pkcrop.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pkdiff.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pkdsm2shadow.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pkdumpimg.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pkdumpogr.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pkegcs.Po@am__quote@ @@ -551,7 +554,6 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pksieve.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pkstat.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pkstatogr.Po@am__quote@ -@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pkveg2shadow.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/svm.Po@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 _______________________________________________ Pkg-grass-devel mailing list Pkg-grass-devel@lists.alioth.debian.org http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-grass-devel