diff --git a/Licenses/libsvm.txt b/Licenses/libsvm.txt new file mode 100644 index 0000000..e3f0889 --- /dev/null +++ b/Licenses/libsvm.txt @@ -0,0 +1,31 @@ + +Copyright (c) 2000-2019 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/Source/SVMWrapper.cpp b/Source/SVMWrapper.cpp new file mode 100644 index 0000000..5424619 --- /dev/null +++ b/Source/SVMWrapper.cpp @@ -0,0 +1,192 @@ +// SVMWrapper.cpp : Defines the exported functions for the DLL application. +// + +#include "stdafx.h" + +#include +#include "svm.h" + +using namespace std; + +extern "C" void __cdecl print_log(const char *str) { + return; +} + +struct SVMWrapper { + struct svm_model *model; + int cols; + double *means; + double *sdevs; + struct svm_parameter param; + struct svm_problem *problem; +}; + +typedef struct svm_node svmnode; + +svmnode **make_svmnode_array(float features[], int rows, int cols, double *means, double *sdevs, int rescale) { + svmnode **x = new svmnode*[rows]; + svmnode *nodes = new svmnode[rows * (cols + 1)]; + + double row_factor = 1.0 / (double)rows; + if (rescale) { + for (int c = 0; c < cols; c++) { + means[c] = 0; + sdevs[c] = 0; + } + for (int r = 0; r < rows; r++) { + for (int c = 0; c < cols; c++) { + means[c] += (double)features[r * cols + c]; + } + } + for (int c = 0; c < cols; c++) { + means[c] *= row_factor; + } + for (int r = 0; r < rows; r++) { + for (int c = 0; c < cols; c++) { + double sq_diff = ((double)features[r * cols + c] - means[c]); + sq_diff *= sq_diff; + sdevs[c] += sq_diff; + } + } + for (int c = 0; c < cols; c++) { + sdevs[c] *= row_factor; + sdevs[c] = sqrt(sdevs[c]); + sdevs[c] = 1.0 / sdevs[c]; + } + } + + for (int r = 0; r < rows; r++) { + x[r] = &(nodes[r * (cols + 1)]); + int idx = 0; + for (int c = 0; c < cols; c++) { + x[r][c].index = c + 1; + x[r][c].value = ((double)features[r * cols + c] - means[c]) * sdevs[c]; + } + x[r][cols].index = -1; + } + return x; +} + +void make_problem(struct svm_problem *problem, float features[], float labels[], int rows, int cols, double *means, double *sdevs) { + problem->l = rows; + problem->y = new double[rows]; + for (int i = 0; i < rows; i++) + problem->y[i] = (double)labels[i]; + problem->x = make_svmnode_array(features, rows, cols, means, sdevs, 1); +} + +void destroy_svmnode_array(svmnode **x) { + delete x[0]; + delete x; +} + +void destroy_problem(struct svm_problem *problem) { + delete problem->y; + destroy_svmnode_array(problem->x); +} + +extern "C" { + __declspec(dllexport) void *__cdecl trainModel(float features[], float labels[], int rows, int cols) { + if (!(rows > 0 && cols > 0)) + return NULL; + + svm_set_print_string_function(print_log); + + struct SVMWrapper *sw = new SVMWrapper(); + sw->cols = cols; + sw->means = new double[cols]; + sw->sdevs = new double[cols]; + + sw->param.svm_type = C_SVC; + sw->param.kernel_type = RBF; + sw->param.degree = 3; + sw->param.gamma = 1.0 / (double)cols; + sw->param.coef0 = 0; + sw->param.nu = 0.5; + sw->param.cache_size = 100; + sw->param.C = 4; + sw->param.eps = 0.001; + sw->param.p = 0.1; + sw->param.shrinking = 1; + sw->param.probability = 0; + sw->param.nr_weight = 0; + sw->param.weight_label = NULL; + sw->param.weight = NULL; + + sw->problem = new struct svm_problem(); + make_problem(sw->problem, features, labels, rows, cols, sw->means, sw->sdevs); + sw->model = svm_train(sw->problem, &sw->param); + + return sw; + } + + __declspec(dllexport) void __cdecl predict(void *ptr, float features[], float predictions[], int rows) { + struct SVMWrapper *sw = (struct SVMWrapper *)ptr; + int cols = sw->cols; + + svmnode **x = make_svmnode_array(features, rows, cols, sw->means, sw->sdevs, 0); + for (int r = 0; r < rows; r++) + predictions[r] = (float)svm_predict(sw->model, x[r]); + destroy_svmnode_array(x); + + return; + } + + __declspec(dllexport) void *__cdecl loadModel(char *filename, int cols, double means[], double sdevs[]) { + struct SVMWrapper *sw = new SVMWrapper(); + sw->cols = cols; + sw->model = svm_load_model(filename); + sw->means = new double[cols]; + sw->sdevs = new double[cols]; + for (int c = 0; c < cols; c++) { + sw->means[c] = means[c]; + sw->sdevs[c] = sdevs[c]; + } + return sw; + } + + __declspec(dllexport) void *__cdecl loadModelString(char *modelString, int cols, double means[], double sdevs[]) { + struct SVMWrapper *sw = new SVMWrapper(); + sw->cols = cols; + sw->model = svm_load_model_string(modelString); + sw->means = new double[cols]; + sw->sdevs = new double[cols]; + for (int c = 0; c < cols; c++) { + sw->means[c] = means[c]; + sw->sdevs[c] = sdevs[c]; + } + sw->problem = NULL; + return sw; + } + + __declspec(dllexport) void __cdecl getScales(void *ptr, double means[], double sdevs[]) { + struct SVMWrapper *sw = (struct SVMWrapper *)ptr; + int cols = sw->cols; + for (int c = 0; c < cols; c++) { + means[c] = sw->means[c]; + sdevs[c] = sw->sdevs[c]; + } + } + + __declspec(dllexport) char *__cdecl getModel(void *ptr) { + struct SVMWrapper *sw = (struct SVMWrapper *)ptr; + return svm_save_model_string(sw->model); + } + + __declspec(dllexport) int __cdecl saveModel(void *ptr, char *filename) { + struct SVMWrapper *sw = (struct SVMWrapper *)ptr; + return svm_save_model(filename, sw->model); + } + + __declspec(dllexport) void __cdecl destroyModel(void *ptr) { + struct SVMWrapper *sw = (struct SVMWrapper *)ptr; + delete sw->means; + delete sw->sdevs; + svm_free_and_destroy_model(&sw->model); + if (sw->problem != NULL) { + destroy_problem(sw->problem); + delete sw->problem; + } + delete sw; + } +} \ No newline at end of file diff --git a/Source/svm.cpp b/Source/svm.cpp new file mode 100644 index 0000000..0271fca --- /dev/null +++ b/Source/svm.cpp @@ -0,0 +1,3583 @@ +#include "stdafx.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "svm.h" +int libsvm_version = LIBSVM_VERSION; +typedef float Qfloat; +typedef signed char schar; +#ifndef min +template static inline T min(T x,T y) { return (x static inline T max(T x,T y) { return (x>y)?x:y; } +#endif +template static inline void swap(T& x, T& y) { T t=x; x=y; y=t; } +template 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;iindex != -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 2*active_size*(l-active_size)) + { + for(i=active_size;iget_Q(i,active_size); + for(j=0;jget_Q(i,l); + double alpha_i = alpha[i]; + for(j=active_size;jl = 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;iINT_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= max_iter) + { + if(active_size < l) + { + // reconstruct the whole gradient to calculate objective value + reconstruct_gradient(); + active_size = l; + info("*"); + } + //fprintf(stderr,"\nWARNING: reaching max number of iterations\n"); + } + + // calculate rho + + si->rho = calculate_rho(); + + // calculate objective value + { + double v = 0; + int i; + for(i=0;iobj = v/2; + } + + // put back the solution + { + for(int i=0;iupper_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= 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= 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 || Gmin_idx == -1) + 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= 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 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;i0) + 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= 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= 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 || Gmin_idx == -1) + 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 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 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 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*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*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*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*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*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*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;jl; + double *minus_ones = new double[l]; + schar *y = new schar[l]; + + int i; + + for(i=0;iy[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;il)); + + for(i=0;il; + double nu = param->nu; + + schar *y = new schar[l]; + + for(i=0;iy[i]>0) + y[i] = +1; + else + y[i] = -1; + + double sum_pos = nu*l/2; + double sum_neg = nu*l/2; + + for(i=0;ieps, si, param->shrinking); + double r = si->r; + + info("C = %f\n",1/r); + + for(i=0;irho /= 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;il) + alpha[n] = param->nu * prob->l - n; + for(i=n+1;ieps, 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;ip - 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;iC*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;iy[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;il); + 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;il;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 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;i0) 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= 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)= min_step) + { + newA = A + stepsize * dA; + newB = B + stepsize * dB; + + // New function value + newf = 0.0; + for (i=0;i= 0) + newf += t[i]*fApB + log(1+exp(-fApB)); + else + newf += (t[i] - 1)*fApB +log(1+exp(fApB)); + } + // Check sufficient decrease + if (newf=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;tmax_error) + max_error=error; + } + if (max_error=max_iter) + info("Exceeds max_iter in multiclass_prob\n"); + for(t=0;tl); + double *dec_values = Malloc(double,prob->l); + + // random shuffle + for(i=0;il;i++) perm[i]=i; + for(i=0;il;i++) + { + int j = i+rand()%(prob->l-i); + swap(perm[i],perm[j]); + } + for(i=0;il/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;jx[perm[j]]; + subprob.y[k] = prob->y[perm[j]]; + ++k; + } + for(j=end;jl;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;j0) + p_count++; + else + n_count++; + + if(p_count==0 && n_count==0) + for(j=begin;j 0 && n_count == 0) + for(j=begin;j 0) + for(j=begin;jx[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;il;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;il;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;iy[i]; + int j; + for(j=0;jparam = *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;il;i++) + if(fabs(f.alpha[i]) > 0) ++nSV; + model->l = nSV; + model->SV = Malloc(svm_node *,nSV); + model->sv_coef[0] = Malloc(double,nSV); + model->sv_indices = Malloc(int,nSV); + int j = 0; + for(i=0;il;i++) + if(fabs(f.alpha[i]) > 0) + { + model->SV[j] = prob->x[i]; + model->sv_coef[0][j] = f.alpha[i]; + model->sv_indices[j] = i+1; + ++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;ix[perm[i]]; + + // calculate weighted C + + double *weighted_C = Malloc(double, nr_class); + for(i=0;iC; + for(i=0;inr_weight;i++) + { + int j; + for(j=0;jweight_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;iprobability) + { + probA=Malloc(double,nr_class*(nr_class-1)/2); + probB=Malloc(double,nr_class*(nr_class-1)/2); + } + + int p = 0; + for(i=0;iprobability) + 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 0) + nonzero[si+k] = true; + for(k=0;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;ilabel[i] = label[i]; + + model->rho = Malloc(double,nr_class*(nr_class-1)/2); + for(i=0;irho[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;iprobA[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;inSV[i] = nSV; + nz_count[i] = nSV; + } + + info("Total nSV = %d\n",total_sv); + + model->l = total_sv; + model->SV = Malloc(svm_node *,total_sv); + model->sv_indices = Malloc(int,total_sv); + p = 0; + for(i=0;iSV[p] = x[i]; + model->sv_indices[p++] = perm[i] + 1; + } + + int *nz_start = Malloc(int,nr_class); + nz_start[0] = 0; + for(i=1;isv_coef = Malloc(double *,nr_class-1); + for(i=0;isv_coef[i] = Malloc(double,total_sv); + + p = 0; + for(i=0;isv_coef[j-1][q++] = f[p].alpha[k]; + q = nz_start[j]; + for(k=0;ksv_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;il; + int *perm = Malloc(int,l); + int nr_class; + if (nr_fold > l) + { + nr_fold = l; + //fprintf(stderr,"WARNING: # folds > # data. Will use # folds = # data instead (i.e., leave-one-out cross validation)\n"); + } + fold_start = Malloc(int,nr_fold+1); + // 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;ix[perm[j]]; + subprob.y[k] = prob->y[perm[j]]; + ++k; + } + for(j=end;jx[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;jx[perm[j]],prob_estimates); + free(prob_estimates); + } + else + for(j=begin;jx[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;inr_class;i++) + label[i] = model->label[i]; +} + +void svm_get_sv_indices(const svm_model *model, int* indices) +{ + if (model->sv_indices != NULL) + for(int i=0;il;i++) + indices[i] = model->sv_indices[i]; +} + +int svm_get_nr_sv(const svm_model *model) +{ + return model->l; +} + +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;il;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;iSV[i],model->param); + + int *start = Malloc(int,nr_class); + start[0] = 0; + for(i=1;inSV[i-1]; + + int *vote = Malloc(int,nr_class); + for(i=0;inSV[i]; + int cj = model->nSV[j]; + + int k; + double *coef1 = model->sv_coef[j-1]; + double *coef2 = model->sv_coef[i]; + for(k=0;krho[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 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;iprobA[k],model->probB[k]),min_prob),1-min_prob); + pairwise_prob[j][i]=1-pairwise_prob[i][j]; + k++; + } + if (nr_class == 2) + { + prob_estimates[0] = pairwise_prob[0][1]; + prob_estimates[1] = pairwise_prob[1][0]; + } + else + multiclass_probability(nr_class,pairwise_prob,prob_estimates); + + int prob_max_idx = 0; + for(i=1;i prob_estimates[prob_max_idx]) + prob_max_idx = i; + for(i=0;ilabel[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 = setlocale(LC_ALL, NULL); + if (old_locale) { + old_locale = _strdup(old_locale); + } + 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 %.17g\n", param.gamma); + + if (param.kernel_type == POLY || param.kernel_type == SIGMOID) + fprintf(fp, "coef0 %.17g\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, " %.17g", 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, " %.17g", 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, " %.17g", 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, "%.17g ", 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; +} + +int vasprintf(char **strp, const char *format, va_list ap) +{ + int len = _vscprintf(format, ap); + if (len == -1) + return -1; + char *str = (char*)malloc((size_t)len + 1); + if (!str) + return -1; + int retval = vsnprintf(str, len + 1, format, ap); + if (retval == -1) { + free(str); + return -1; + } + *strp = str; + return retval; +} + +int asprintf(char **strp, const char *format, ...) +{ + va_list ap; + va_start(ap, format); + int retval = vasprintf(strp, format, ap); + va_end(ap); + return retval; +} + +char *svm_save_model_string(const svm_model *model) +{ + char *old_locale = setlocale(LC_ALL, NULL); + if (old_locale) { + old_locale = _strdup(old_locale); + } + setlocale(LC_ALL, "C"); + + std::vector strings = std::vector(); + char *tmp; + + const svm_parameter& param = model->param; + + asprintf(&tmp, "svm_type %s\n", svm_type_table[param.svm_type]); + strings.push_back(tmp); + asprintf(&tmp, "kernel_type %s\n", kernel_type_table[param.kernel_type]); + strings.push_back(tmp); + + if (param.kernel_type == POLY) { + asprintf(&tmp, "degree %d\n", param.degree); + strings.push_back(tmp); + } + + if (param.kernel_type == POLY || param.kernel_type == RBF || param.kernel_type == SIGMOID) { + asprintf(&tmp, "gamma %.17g\n", param.gamma); + strings.push_back(tmp); + } + + if (param.kernel_type == POLY || param.kernel_type == SIGMOID) { + asprintf(&tmp, "coef0 %.17g\n", param.coef0); + strings.push_back(tmp); + } + + int nr_class = model->nr_class; + int l = model->l; + asprintf(&tmp, "nr_class %d\n", nr_class); + strings.push_back(tmp); + asprintf(&tmp, "total_sv %d\n", l); + strings.push_back(tmp); + + { + asprintf(&tmp, "rho"); + strings.push_back(tmp); + for (int i = 0; i < nr_class*(nr_class - 1) / 2; i++) { + asprintf(&tmp, " %.17g", model->rho[i]); + strings.push_back(tmp); + } + asprintf(&tmp, "\n"); + strings.push_back(tmp); + } + + if (model->label) + { + asprintf(&tmp, "label"); + strings.push_back(tmp); + for (int i = 0; i < nr_class; i++) { + asprintf(&tmp, " %d", model->label[i]); + strings.push_back(tmp); + } + asprintf(&tmp, "\n"); + strings.push_back(tmp); + } + + if (model->probA) // regression has probA only + { + asprintf(&tmp, "probA"); + strings.push_back(tmp); + for (int i = 0; i < nr_class*(nr_class - 1) / 2; i++) { + asprintf(&tmp, " %.17g", model->probA[i]); + strings.push_back(tmp); + } + asprintf(&tmp, "\n"); + strings.push_back(tmp); + } + if (model->probB) + { + asprintf(&tmp, "probB"); + strings.push_back(tmp); + for (int i = 0; i < nr_class*(nr_class - 1) / 2; i++) { + asprintf(&tmp, " %.17g", model->probB[i]); + strings.push_back(tmp); + } + asprintf(&tmp, "\n"); + strings.push_back(tmp); + } + + if (model->nSV) + { + asprintf(&tmp, "nr_sv"); + strings.push_back(tmp); + for (int i = 0; i < nr_class; i++) { + asprintf(&tmp, " %d", model->nSV[i]); + strings.push_back(tmp); + } + asprintf(&tmp, "\n"); + strings.push_back(tmp); + } + + asprintf(&tmp, "SV\n"); + strings.push_back(tmp); + 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++) { + asprintf(&tmp, "%.17g ", sv_coef[j][i]); + strings.push_back(tmp); + } + + const svm_node *p = SV[i]; + + if (param.kernel_type == PRECOMPUTED) { + asprintf(&tmp, "0:%d ", (int)(p->value)); + strings.push_back(tmp); + } + else + while (p->index != -1) + { + asprintf(&tmp, "%d:%.8g ", p->index, p->value); + strings.push_back(tmp); + p++; + } + asprintf(&tmp, "\n"); + strings.push_back(tmp); + } + + size_t length = 0; + for (int i = 0; i < strings.size(); i++) { + length += strlen(strings[i]); + } + char *save = (char *)CoTaskMemAlloc(length + 1); + save[length] = 0; + char *save_ptr = save; + for (int i = 0; i < strings.size(); i++) { + size_t len = strlen(strings[i]); + memcpy(save_ptr, strings[i], len); + save_ptr += len; + free(strings[i]); + } + + setlocale(LC_ALL, old_locale); + free(old_locale); + + return save; +} + +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; +} + +// +// FSCANF helps to handle fscanf failures. +// Its do-while block avoids the ambiguity when +// if (...) +// FSCANF(); +// is used +// +#define FSCANF(_stream, _format, _var) do{ if (fscanf(_stream, _format, _var) != 1) return false; }while(0) +bool read_model_header(FILE *fp, svm_model* model) +{ + svm_parameter& param = model->param; + // parameters for training only won't be assigned, but arrays are assigned as NULL for safety + param.nr_weight = 0; + param.weight_label = NULL; + param.weight = 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"); + return false; + } + } + 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"); + return false; + } + } + 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;irho[i]); + } + else if(strcmp(cmd,"label")==0) + { + int n = model->nr_class; + model->label = Malloc(int,n); + for(int i=0;ilabel[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;iprobA[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;iprobB[i]); + } + else if(strcmp(cmd,"nr_sv")==0) + { + int n = model->nr_class; + model->nSV = Malloc(int,n); + for(int i=0;inSV[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); + return false; + } + } + + return true; + +} + +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 = setlocale(LC_ALL, NULL); + if (old_locale) { + old_locale = _strdup(old_locale); + } + setlocale(LC_ALL, "C"); + + // read parameters + + svm_model *model = Malloc(svm_model,1); + model->rho = NULL; + model->probA = NULL; + model->probB = NULL; + model->sv_indices = NULL; + model->label = NULL; + model->nSV = NULL; + + // read header + if (!read_model_header(fp, model)) + { + //fprintf(stderr, "ERROR: fscanf failed to read model\n"); + 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;isv_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;iSV[i] = &x_space[j]; + + p = strtok(line, " \t"); + model->sv_coef[0][i] = strtod(p,&endptr); + for(int k=1;ksv_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; +} + +#define SCANF(_stream, _format, _var) do{ int scan_pos; if (sscanf(*_stream, _format "%n", _var, &scan_pos) != 1) return false; *_stream += scan_pos; }while(0) +bool read_model_header_string(char **model_string, svm_model* model) +{ + svm_parameter& param = model->param; + // parameters for training only won't be assigned, but arrays are assigned as NULL for safety + param.nr_weight = 0; + param.weight_label = NULL; + param.weight = NULL; + + char cmd[81]; + while (1) + { + SCANF(model_string, "%80s", cmd); + + if (strcmp(cmd, "svm_type") == 0) + { + SCANF(model_string, "%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"); + return false; + } + } + else if (strcmp(cmd, "kernel_type") == 0) + { + SCANF(model_string, "%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"); + return false; + } + } + else if (strcmp(cmd, "degree") == 0) + SCANF(model_string, "%d", ¶m.degree); + else if (strcmp(cmd, "gamma") == 0) + SCANF(model_string, "%lf", ¶m.gamma); + else if (strcmp(cmd, "coef0") == 0) + SCANF(model_string, "%lf", ¶m.coef0); + else if (strcmp(cmd, "nr_class") == 0) + SCANF(model_string, "%d", &model->nr_class); + else if (strcmp(cmd, "total_sv") == 0) + SCANF(model_string, "%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++) + SCANF(model_string, "%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++) + SCANF(model_string, "%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++) + SCANF(model_string, "%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++) + SCANF(model_string, "%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++) + SCANF(model_string, "%d", &model->nSV[i]); + } + else if (strcmp(cmd, "SV") == 0) + { + while (1) + { + int c = *((*model_string)++); + if (c == EOF || c == '\n') break; + } + break; + } + else + { + //fprintf(stderr,"unknown text in model file: [%s]\n",cmd); + return false; + } + } + + return true; + +} + +struct svm_model *svm_load_model_string(char *model_string) { + char *old_locale = setlocale(LC_ALL, NULL); + if (old_locale) { + old_locale = _strdup(old_locale); + } + setlocale(LC_ALL, "C"); + + // read parameters + + svm_model *model = Malloc(svm_model, 1); + model->rho = NULL; + model->probA = NULL; + model->probB = NULL; + model->sv_indices = NULL; + model->label = NULL; + model->nSV = NULL; + + // read header + if (!read_model_header_string(&model_string, model)) + { + //fprintf(stderr, "ERROR: fscanf failed to read model\n"); + 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; + char *pos = model_string; + + max_line_len = 1024; + line = Malloc(char, strlen(model_string) + 2); + memcpy(line, model_string, strlen(model_string) + 1); + char *p, *endptr, *idx, *val; + + //while (readline_string(&model_string) != 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); + size_t mlen = strlen(model_string); + int mend = 0; + for (int i = 0; i < mlen && model_string[i] != '\n' && model_string[i] != '\0'; i++) + mend++; + memcpy(line, model_string, mend + 1); + line[mend + 2] = 0; + model_string += mend + 1; + + 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); + + 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;inr_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->sv_indices); + model_ptr->sv_indices = 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((kernel_type == POLY || kernel_type == RBF || kernel_type == SIGMOID) && + param->gamma < 0) + return "gamma < 0"; + + if(kernel_type == POLY && 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;iy[i]; + int j; + for(j=0;jnu*(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/Source/svm.h b/Source/svm.h new file mode 100644 index 0000000..6d184b4 --- /dev/null +++ b/Source/svm.h @@ -0,0 +1,106 @@ +#ifndef _LIBSVM_H +#define _LIBSVM_H + +#define LIBSVM_VERSION 324 + +#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; + int *sv_indices; /* sv_indices[0,...,nSV-1] are values in [1,...,num_traning_data] to indicate SVs in the training set */ + + /* 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); +char *svm_save_model_string(const svm_model *model); +struct svm_model *svm_load_model(const char *model_file_name); +struct svm_model *svm_load_model_string(char *model_string); + +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); +void svm_get_sv_indices(const struct svm_model *model, int *sv_indices); +int svm_get_nr_sv(const struct svm_model *model); +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/Unity/OpenSee.cs b/Unity/OpenSee.cs index fd70b3c..70c4e70 100644 --- a/Unity/OpenSee.cs +++ b/Unity/OpenSee.cs @@ -16,7 +16,7 @@ public class OpenSee : MonoBehaviour { public int listenPort = 11573; private const int nPoints = 66; - private const int packetFrameSize = 8 + 2 * 4 + 2 * 4 + 1 + 4 + 3 * 4 + 3 * 4 + 4 * 4 + 4 * 68 + 4 * 2 * 68 + 4 * 3 * 70; + private const int packetFrameSize = 8 + 4 + 2 * 4 + 2 * 4 + 1 + 4 + 3 * 4 + 3 * 4 + 4 * 4 + 4 * 68 + 4 * 2 * 68 + 4 * 3 * 70; [Header("Tracking data")] [Tooltip("This is an informational property that tells you how many packets have been received")] @@ -28,7 +28,7 @@ public class OpenSee : MonoBehaviour { public class OpenSeeData { [Tooltip("The time this tracking data was captured at.")] public double time; - [Tooltip("This is the id of the tracked face. When tracking multiple faces, they might get reordered due to faces coming and going, but as long as tracking is not lost on a face, its id should stay the same.")] + [Tooltip("This is the id of the tracked face. When tracking multiple faces, they might get reordered due to faces coming and going, but as long as tracking is not lost on a face, its id should stay the same. Face ids depend only on the order of first detection and locations of the faces.")] public int id; [Tooltip("This field gives the resolution of the camera or video being tracked.")] public Vector2 cameraResolution; diff --git a/Unity/OpenSee.cs.meta b/Unity/OpenSee.cs.meta new file mode 100644 index 0000000..f870d69 --- /dev/null +++ b/Unity/OpenSee.cs.meta @@ -0,0 +1,11 @@ +fileFormatVersion: 2 +guid: 7b5af009eeb86f9469d2be8970e5a0f4 +MonoImporter: + externalObjects: {} + serializedVersion: 2 + defaultReferences: [] + executionOrder: 0 + icon: {instanceID: 0} + userData: + assetBundleName: + assetBundleVariant: diff --git a/Unity/OpenSeeExpression.cs b/Unity/OpenSeeExpression.cs new file mode 100644 index 0000000..785d71f --- /dev/null +++ b/Unity/OpenSeeExpression.cs @@ -0,0 +1,406 @@ +using System; +using System.Collections.Generic; +using System.IO; +using System.IO.Compression; +using System.Runtime.Serialization; +using System.Runtime.Serialization.Formatters.Binary; +using UnityEngine; + +namespace OpenSee { + +public class OpenSeeExpression : MonoBehaviour +{ + [Header("Settings")] + [Tooltip("This is the source of face tracking data. If it is not set, any OpenSee component on the current object is used.")] + public OpenSee openSee; + [Tooltip("This specifies which expression calibration data should be collected for.")] + public string calibrationExpression = "neutral"; + [Tooltip("This specifies the id of the face for which data should be collected and predictions should be made. Face ids depend only on the order of first detection and locations of the faces.")] + public int faceId = 0; + [Tooltip("This keeps low numbers of misdetections from changing the expression, but introduces a bit of lag to expression changes. For example, a value of 1 means that a new expression has to be detected twice in a row, before the predicted expression changes.")] + public int expressionStabilizer = 1; + [Tooltip("Setting this to a value above 1 will record only every recordingSkip-th frame when capture calibration data. This can be useful with cameras that have a high framerate leading to only a short time for capturing varied expression information. Usually at least 20-30 seconds are required to capture a good variety of angles and expression variation, so this should be set accordingly when using the recording flag. At 1, no frames are skipped.")] + public int recordingSkip = 2; + [Tooltip("This is the filename used for loading and saving expression data using the load and save flags.")] + public string filename = ""; + [Header("Toggles")] + [Tooltip("When enabled, calibration data will be collected from the given OpenSee component.")] + public bool recording = false; + [Tooltip("When enabled, the calibration data collected for the current expression is cleared and this flag is set back to false.")] + public bool clear = false; + [Tooltip("When enabled, a new prediction model is trained and this flag is set back to false. Only expressions for which percentRecorded is 100% will be included in training.")] + public bool train = false; + [Tooltip("When enabled and a model is loaded, the current expression will be predicted every frame.")] + public bool predict = false; + [Tooltip("When enabled, the expression data from the specified filename will be loaded and this flag is set back to false.")] + public bool load = false; + [Tooltip("When enabled, the expression data will be saved to the specified filename and this flag is set back to false.")] + public bool save = false; + [Header("Information")] + [Tooltip("This is the of expressions for which calibration data was collected. The maximum number is 25.")] + public int expressionNumber = 0; + [Tooltip("This is the percentage of necessary training data collected for the current expression.")] + public float percentRecorded = 0.0f; + [Tooltip("This shows whether a model is ready to predict expressions with.")] + public bool modelReady = false; + [Tooltip("This is the time the current expression was detected.")] + public float expressionTime = 0f; + [Tooltip("This is currently detected expression.")] + public string expression = ""; + [Header("Training statistics")] + [Tooltip("This shows the accuracy result of the last training run.")] + public float accuracy = 0.0f; + [Tooltip("This is the confusion matrix for the last training run's validation set.")] + public int[,] confusionMatrix = null; + [Tooltip("This is a pretty printed string representation of the confusion matrix.")] + public string confusionMatrixString = ""; + [Tooltip("This shows any accuracy warnings that resulted from the last training run.")] + public string[] warnings = null; + + [Serializable] + private class OpenSeeExpressionRepresentation { + private Dictionary> expressions; + private byte[] modelBytes = null; + private string[] classLabels = null; + + static public void LoadSerialized(byte[] modelBytes, out Dictionary> expressions, out SVMModel model, out string[] classLabels) { + IFormatter formatter = new BinaryFormatter(); + MemoryStream memoryStream = new MemoryStream(modelBytes); + memoryStream.Position = 0; + OpenSeeExpressionRepresentation oser; + using (GZipStream gzipStream = new GZipStream(memoryStream, CompressionMode.Decompress)) { + oser = formatter.Deserialize(gzipStream) as OpenSeeExpressionRepresentation; + } + expressions = oser.expressions; + model = new SVMModel(oser.modelBytes); + classLabels = oser.classLabels; + } + + static public byte[] ToSerialized(Dictionary> expressions, SVMModel model, string[] classLabels) { + OpenSeeExpressionRepresentation oser = new OpenSeeExpressionRepresentation(); + oser.expressions = expressions; + oser.modelBytes = model.SaveModel(); + oser.classLabels = classLabels; + + IFormatter formatter = new BinaryFormatter(); + MemoryStream memoryStream = new MemoryStream(); + using (GZipStream gzipStream = new GZipStream(memoryStream, CompressionMode.Compress)) { + formatter.Serialize(gzipStream, oser); + gzipStream.Flush(); + } + return memoryStream.ToArray(); + } + } + + private Dictionary> expressions; + private SVMModel model = null; + private string[] classLabels = null; + private int maxSamples = 400; // Allows 25 expressions without exceeding 10000 rows in total for SVM training, which is a commonly seen upper limit. + // rightEyeOpen, leftEyeOpen, translation, rawQuaternion, rawEuler, confidence, points/(width, height), points3D + private int cols = 1 + 1 + 3 + 4 + 3 + /*66 + 2 * 66 +*/ 3 * 66; + private double lastCapture = 0.0; + private float warningThreshold = 5; + private int currentPrediction = -1; + private int lastPrediction = -1; + private int lastPredictionCount = 0; + private int frameCount = 0; + + private void ResetInfo() { + recording = false; + clear = false; + train = false; + predict = false; + modelReady = false; + accuracy = 0.0f; + confusionMatrix = null; + warnings = null; + percentRecorded = 0.0f; + expression = ""; + confusionMatrixString = ""; + lastPrediction = -1; + currentPrediction = -1; + lastPredictionCount = 0; + expressionTime = 0f; + } + + void Start() + { + if (openSee == null) + openSee = GetComponent(); + if (openSee == null) { + ResetInfo(); + return; + } + ResetInfo(); + expressions = new Dictionary>(); + model = new SVMModel(); + } + + private float[] GetData(OpenSee.OpenSeeData t) { + float[] data = new float[cols]; + data[0] = t.rightEyeOpen; + data[1] = t.leftEyeOpen; + data[2] = t.translation.x; + data[3] = t.translation.y; + data[4] = t.translation.z; + data[5] = t.rawQuaternion.x; + data[6] = t.rawQuaternion.y; + data[7] = t.rawQuaternion.z; + data[8] = t.rawQuaternion.w; + data[9] = t.rawEuler.x; + data[10] = t.rawEuler.y; + data[11] = t.rawEuler.z; + /*for (int i = 0; i < 66; i++) + data[12 + i] = t.confidence[i]; + for (int i = 0; i < 66; i++) { + data[12 + 66 + 2 * i] = t.points[i].x / t.cameraResolution.x; + data[12 + 66 + 2 * i + 1] = t.points[i].y / t.cameraResolution.y; + }*/ + for (int i = 0; i < 66; i++) { + data[12 + /*66 + 2 * 66 +*/ 3 * i] = t.points3D[i].x; + data[12 + /*66 + 2 * 66 +*/ 3 * i + 1] = t.points3D[i].y; + data[12 + /*66 + 2 * 66 +*/ 3 * i + 2] = t.points3D[i].z; + } + return data; + } + + public bool TrainModel() { + train = false; + if (openSee == null) { + ResetInfo(); + return false; + } + List keys = new List(); + foreach (string key in expressions.Keys) { + List list = expressions[key]; + if (list.Count == maxSamples) + keys.Add(key); + else + return false; + } + int classes = keys.Count; + if (classes < 2 || classes > 25) + return false; + keys.Sort(); + classLabels = keys.ToArray(); + int train_split = maxSamples * 3 / 4; + int test_split = maxSamples - train_split; + int rows_train = classes * train_split; + int rows_test = classes * test_split; + float[] X_train = new float[rows_train * cols]; + float[] X_test = new float[rows_test * cols]; + float[] y_train = new float[rows_train]; + float[] y_test = new float[rows_test]; + int i_train = 0; + int i_test = 0; + System.Random rnd = new System.Random(); + for (int i = 0; i < classes; i++) { + List list = expressions[keys[i]]; + for (int j = maxSamples; j > 1;) { + j--; + int k = rnd.Next(j + 1); + float[] tmp = list[k]; + list[k] = list[j]; + list[j]= tmp; + } + for (int j = 0; j < train_split; j++) { + for (int k = 0; k < cols; k++) { + X_train[i_train * cols + k] = list[j][k]; + } + y_train[i_train] = i; + i_train++; + } + for (int j = train_split; j < train_split + test_split; j++) { + for (int k = 0; k < cols; k++) + X_test[i_test * cols + k] = list[j][k]; + y_test[i_test] = i; + i_test++; + } + } + model.TrainModel(X_train, y_train, rows_train, cols); + confusionMatrix = model.ConfusionMatrix(X_test, y_test, rows_test, out accuracy); + confusionMatrixString = SVMModel.FormatMatrix(confusionMatrix, classLabels); + List accuracyWarnings = new List(); + for (int label = 0; label < classes; label++) { + float error = 100f * (1f - ((float)confusionMatrix[label, label] / (float)test_split)); + if (error > warningThreshold) + accuracyWarnings.Add("The expression \"" + classLabels[label] + "\" is misclassified with a chance of " + error.ToString("0.00") + "%."); + } + warnings = accuracyWarnings.ToArray(); + modelReady = true; + return true; + } + + public bool PredictExpression() { + if (openSee == null || !(modelReady && model.Ready())) { + ResetInfo(); + return false; + } + OpenSee.OpenSeeData[] openSeeData = openSee.trackingData; + if (openSeeData == null) + return false; + OpenSee.OpenSeeData t = null; + foreach (OpenSee.OpenSeeData data in openSeeData) + if (data.id == faceId) + t = data; + if (t == null || t.time <= lastCapture) + return false; + float[] faceData = GetData(t); + float[] prediction = model.Predict(faceData, 1); + int predictedExpression = (int)Mathf.Round(prediction[0]); + if (predictedExpression == lastPrediction) + lastPredictionCount++; + else { + lastPrediction = predictedExpression; + lastPredictionCount = 1; + } + if (lastPredictionCount > expressionStabilizer) { + currentPrediction = lastPrediction; + expressionTime = Time.time; + } + if (currentPrediction >= 0 && currentPrediction <= classLabels.Length) + expression = classLabels[currentPrediction]; + else + expression = ""; + return true; + } + + public bool CaptureCalibrationData() { + if (openSee == null) { + ResetInfo(); + return false; + } + if (recordingSkip < 1) + recordingSkip = 1; + if (calibrationExpression == "") + return false; + if (expressions.ContainsKey(calibrationExpression) || expressions.Keys.Count < 25) { + OpenSee.OpenSeeData[] openSeeData = openSee.trackingData; + if (openSeeData == null) + return false; + OpenSee.OpenSeeData t = null; + foreach (OpenSee.OpenSeeData data in openSeeData) + if (data.id == faceId) + t = data; + if (t == null || t.time <= lastCapture) + return false; + if (!expressions.ContainsKey(calibrationExpression)) { + expressions.Add(calibrationExpression, new List()); + } + List list = expressions[calibrationExpression]; + if (list.Count >= maxSamples) { + percentRecorded = 100f; + recording = false; + return false; + } + if (frameCount % recordingSkip != 0) + return true; + list.Add(GetData(t)); + percentRecorded = 100f * (float)list.Count/(float)maxSamples; + return true; + } else { + recording = false; + calibrationExpression = ""; + return false; + } + } + + public string[] GetExpressions() { + List keys = new List(); + foreach (string key in expressions.Keys) + keys.Add(key); + keys.Sort(); + return keys.ToArray(); + } + + public void ClearExpression() { + if (expressions.ContainsKey(calibrationExpression)) + expressions.Remove(calibrationExpression); + clear = false; + } + + public void LoadFromBytes(byte[] data) { + load = false; + if (openSee == null) { + ResetInfo(); + return; + } + OpenSeeExpressionRepresentation.LoadSerialized(data, out expressions, out model, out classLabels); + if (model.Ready()) + modelReady = true; + } + + public byte[] SaveToBytes() { + save = false; + if (openSee == null) { + ResetInfo(); + return null; + } + return OpenSeeExpressionRepresentation.ToSerialized(expressions, model, classLabels); + } + + public void LoadFromFile() { + load = false; + if (openSee == null) { + ResetInfo(); + return; + } + byte[] data = File.ReadAllBytes(filename); + LoadFromBytes(data); + } + + public void SaveToFile() { + save = false; + if (openSee == null) { + ResetInfo(); + return; + } + byte[] data = SaveToBytes(); + File.WriteAllBytes(filename, data); + } + + void Update() + { + frameCount++; + if (openSee == null || expressions == null) { + ResetInfo(); + return; + } + if (clear) + ClearExpression(); + if (recording && calibrationExpression != "") { + CaptureCalibrationData(); + } + if (train) + TrainModel(); + expressionNumber = expressions.Keys.Count; + if (load) + LoadFromFile(); + if (save) + SaveToFile(); + if (!predict) + expression = ""; + else if (modelReady) + PredictExpression(); + else + predict = false; + + if (expressions.ContainsKey(calibrationExpression)) + percentRecorded = 100f * (float)expressions[calibrationExpression].Count/(float)maxSamples; + else + percentRecorded = 0f; + + OpenSee.OpenSeeData[] openSeeData = openSee.trackingData; + if (openSeeData == null) + return; + OpenSee.OpenSeeData t = null; + foreach (OpenSee.OpenSeeData data in openSeeData) + if (data.id == faceId) + t = data; + if (t == null || t.time <= lastCapture) + return; + lastCapture = t.time; + } +} + +} \ No newline at end of file diff --git a/Unity/OpenSeeExpression.cs.meta b/Unity/OpenSeeExpression.cs.meta new file mode 100644 index 0000000..5b53490 --- /dev/null +++ b/Unity/OpenSeeExpression.cs.meta @@ -0,0 +1,11 @@ +fileFormatVersion: 2 +guid: 9066de4c7793c1e40b212724f72fc064 +MonoImporter: + externalObjects: {} + serializedVersion: 2 + defaultReferences: [] + executionOrder: 0 + icon: {instanceID: 0} + userData: + assetBundleName: + assetBundleVariant: diff --git a/Unity/OpenSeeJobManagement.cs.meta b/Unity/OpenSeeJobManagement.cs.meta new file mode 100644 index 0000000..f522497 --- /dev/null +++ b/Unity/OpenSeeJobManagement.cs.meta @@ -0,0 +1,11 @@ +fileFormatVersion: 2 +guid: 3251308a085a78545a82e7f34ebc2d34 +MonoImporter: + externalObjects: {} + serializedVersion: 2 + defaultReferences: [] + executionOrder: 0 + icon: {instanceID: 0} + userData: + assetBundleName: + assetBundleVariant: diff --git a/Unity/OpenSeeLauncher.cs.meta b/Unity/OpenSeeLauncher.cs.meta new file mode 100644 index 0000000..b707bbd --- /dev/null +++ b/Unity/OpenSeeLauncher.cs.meta @@ -0,0 +1,11 @@ +fileFormatVersion: 2 +guid: beae834477f91834691fb7248fddacf6 +MonoImporter: + externalObjects: {} + serializedVersion: 2 + defaultReferences: [] + executionOrder: 0 + icon: {instanceID: 0} + userData: + assetBundleName: + assetBundleVariant: diff --git a/Unity/OpenSeeShowPoints.cs.meta b/Unity/OpenSeeShowPoints.cs.meta new file mode 100644 index 0000000..cd96b05 --- /dev/null +++ b/Unity/OpenSeeShowPoints.cs.meta @@ -0,0 +1,11 @@ +fileFormatVersion: 2 +guid: 9c0b8bb971ed0e14288631f4274a3807 +MonoImporter: + externalObjects: {} + serializedVersion: 2 + defaultReferences: [] + executionOrder: 0 + icon: {instanceID: 0} + userData: + assetBundleName: + assetBundleVariant: diff --git a/Unity/Plugins.meta b/Unity/Plugins.meta new file mode 100644 index 0000000..9f38e3f --- /dev/null +++ b/Unity/Plugins.meta @@ -0,0 +1,8 @@ +fileFormatVersion: 2 +guid: f1dfc5f1c2e6f294e8191828595882b1 +folderAsset: yes +DefaultImporter: + externalObjects: {} + userData: + assetBundleName: + assetBundleVariant: diff --git a/Unity/Plugins/SVMModel.cs b/Unity/Plugins/SVMModel.cs new file mode 100644 index 0000000..34620fe --- /dev/null +++ b/Unity/Plugins/SVMModel.cs @@ -0,0 +1,248 @@ +using System; +using System.Collections.Generic; +using System.IO; +using System.IO.Compression; +using System.Runtime.InteropServices; +using System.Runtime.Serialization; +using System.Runtime.Serialization.Formatters.Binary; +using System.Text; + +namespace OpenSee { + +public class SVMModel { + #region DllImport + [DllImport("SVMModel", CallingConvention = CallingConvention.Cdecl)] + private static extern System.IntPtr trainModel(float[] features, float[] labels, int rows, int cols); + + [DllImport("SVMModel", CallingConvention = CallingConvention.Cdecl)] + private static extern void predict(System.IntPtr model, float[] features, [Out] float[] predictions, int rows); + + [DllImport("SVMModel", CallingConvention = CallingConvention.Cdecl, CharSet = CharSet.Ansi)] + private static extern System.IntPtr loadModel([MarshalAs(UnmanagedType.LPStr)]string filename, int cols, double[] means, double[] sdevs); + + [DllImport("SVMModel", CallingConvention = CallingConvention.Cdecl, CharSet = CharSet.Ansi)] + private static extern System.IntPtr loadModelString([MarshalAs(UnmanagedType.LPStr)]string modelString, int cols, double[] means, double[] sdevs); + + [DllImport("SVMModel", CallingConvention = CallingConvention.Cdecl)] + private static extern void getScales(System.IntPtr model, [Out] double[] means, [Out] double[] sdevs); + + [DllImport("SVMModel", CallingConvention = CallingConvention.Cdecl, CharSet = CharSet.Ansi)] + private static extern int saveModel(System.IntPtr model, [MarshalAs(UnmanagedType.LPStr)]string filename); + + [DllImport("SVMModel", CallingConvention = CallingConvention.Cdecl, CharSet = CharSet.Ansi)] + private static extern string saveModelString(System.IntPtr model); + + [DllImport("SVMModel", CallingConvention = CallingConvention.Cdecl)] + private static extern void destroyModel(System.IntPtr model); + #endregion + + private bool haveModel = false; + private System.IntPtr model; + private int cols; + private int maxClasses; + + [Serializable] + private class SVMModelRepresentation { + private string modelString; + private int cols; + private int maxClasses; + private double[] means; + private double[] sdevs; + + static public System.IntPtr LoadSerialized(byte[] modelBytes, out int cols, out int maxClasses, bool compress) { + IFormatter formatter = new BinaryFormatter(); + MemoryStream memoryStream = new MemoryStream(modelBytes); + memoryStream.Position = 0; + SVMModelRepresentation smr; + if (compress) { + using (GZipStream gzipStream = new GZipStream(memoryStream, CompressionMode.Decompress)) { + smr = formatter.Deserialize(gzipStream) as SVMModelRepresentation; + } + } else { + smr = formatter.Deserialize(memoryStream) as SVMModelRepresentation; + } + cols = smr.cols; + maxClasses = smr.maxClasses; + System.IntPtr model = loadModelString(smr.modelString, smr.cols, smr.means, smr.sdevs); + return model; + } + + static public byte[] ToSerialized(System.IntPtr model, int cols, int maxClasses, bool compress) { + SVMModelRepresentation smr = new SVMModelRepresentation(); + smr.cols = cols; + smr.maxClasses = maxClasses; + smr.means = new double[cols]; + smr.sdevs = new double[cols]; + getScales(model, smr.means, smr.sdevs); + smr.modelString = saveModelString(model); + + IFormatter formatter = new BinaryFormatter(); + MemoryStream memoryStream = new MemoryStream(); + if (compress) { + using (GZipStream gzipStream = new GZipStream(memoryStream, CompressionMode.Compress)) { + formatter.Serialize(gzipStream, smr); + gzipStream.Flush(); + } + } else { + formatter.Serialize(memoryStream, smr); + } + return memoryStream.ToArray(); + } + } + + public void DestroyModel() { + if (haveModel) { + destroyModel(model); + haveModel = false; + } + } + + private void LoadModel(byte[] modelBytes, bool compress) { + DestroyModel(); + model = SVMModelRepresentation.LoadSerialized(modelBytes, out cols, out maxClasses, compress); + haveModel = true; + } + + public void LoadModel(byte[] modelBytes) { + LoadModel(modelBytes, true); + } + + private byte[] SaveModel(bool compress) { + if (!haveModel) + return null; + return SVMModelRepresentation.ToSerialized(model, cols, maxClasses, compress); + } + + public byte[] SaveModel() { + return SaveModel(true); + } + + // features is a rows*cols long float array with rows written one after another. labels is a rows long float array with the corresponding classes as integral numbers from 0 to num_classes-1, where all classes have to exist. Training on more than 10000 rows is disabled. + public bool TrainModel(float[] features, float[] labels, int rows, int cols) { + DestroyModel(); + if (rows > 10000) + return false; + if (cols < 1 || rows < 1 || features.Length != rows * cols || labels.Length != rows) + return false; + + int max = -1; + HashSet knownClasses = new HashSet(); + for (int i = 0; i < labels.Length; i++) { + labels[i] = (float)Math.Round(labels[i]); + int l = (int)labels[i]; + if (l > max) + max = l; + if (l < 0) + return false; + knownClasses.Add(l); + } + for (int i = 0; i <= max; i++) { + if (!knownClasses.Contains(i)) + return false; + } + + maxClasses = max + 1; + this.cols = cols; + model = trainModel(features, labels, rows, cols); + haveModel = true; + + return true; + } + + // This function takes the same arguments as train, except for cols, which is not necessary. The diagonal are correctly classified entries. The first index is labels, the second index is classifications. + public int[,] ConfusionMatrix(float[] features, float[] labels, int rows, out float accuracy) { + accuracy = 0f; + if (!haveModel) + return null; + if (rows < 1 || features.Length != rows * cols || labels.Length != rows) + return null; + for (int i = 0; i < labels.Length; i++) { + labels[i] = (float)Math.Round(labels[i]); + if (labels[i] < 0 || labels[i] >= maxClasses) + return null; + } + int[,] confusionMatrix = new int[maxClasses, maxClasses]; + float[] predictions = new float[rows]; + predict(model, features, predictions, rows); + for (int i = 0; i < rows; i++) { + predictions[i] = (float)Math.Round(predictions[i]); + if (predictions[i] >= 0 && predictions[i] < maxClasses) { + confusionMatrix[(int)labels[i], (int)predictions[i]]++; + if ((int)labels[i] == (int)predictions[i]) + accuracy++; + } + } + accuracy /= (float)rows; + return confusionMatrix; + } + + static public string FormatMatrix(int[,] matrix) { + return FormatMatrix(matrix, null); + } + + static public string FormatMatrix(int[,] matrix, string[] labels) { + StringBuilder sb = new StringBuilder(); + int rowLength = matrix.GetLength(0); + int colLength = matrix.GetLength(1); + int maxLableLen = 5; + + if (labels != null) { + if (labels.Length != rowLength || labels.Length != colLength) + labels = null; + else { + foreach (string label in labels) + if (label.Length > maxLableLen) + maxLableLen = label.Length; + sb.Append("pred:".PadLeft(maxLableLen)); + sb.Append(" "); + foreach (string label in labels) { + sb.Append(label.PadLeft(maxLableLen)); + sb.Append(" "); + } + sb.Append(Environment.NewLine); + } + } + + for (int r = 0; r < rowLength; r++) { + if (labels != null) { + sb.Append(labels[r].PadLeft(maxLableLen)); + sb.Append(" "); + } + for (int c = 0; c < colLength; c++) { + string n = matrix[r,c].ToString(); + sb.Append(n.PadLeft(maxLableLen)); + sb.Append(" "); + } + sb.Append(Environment.NewLine); + } + return sb.ToString(); + } + + // features is a rows*cols long float array with rows written one after another. + public float[] Predict(float[] features, int rows) { + if (!haveModel) + return null; + if (rows < 1 || features.Length != rows * cols) + return null; + float[] predictions = new float[rows]; + predict(model, features, predictions, rows); + return predictions; + } + + public bool Ready() { + return haveModel; + } + + public SVMModel() { + } + + public SVMModel(byte[] modelBytes) { + LoadModel(modelBytes); + } + + ~SVMModel() { + DestroyModel(); + } +} + +} \ No newline at end of file diff --git a/Unity/Plugins/SVMModel.cs.meta b/Unity/Plugins/SVMModel.cs.meta new file mode 100644 index 0000000..3fa197d --- /dev/null +++ b/Unity/Plugins/SVMModel.cs.meta @@ -0,0 +1,11 @@ +fileFormatVersion: 2 +guid: a95b57bbe200b004dab9b21f0854082e +MonoImporter: + externalObjects: {} + serializedVersion: 2 + defaultReferences: [] + executionOrder: 0 + icon: {instanceID: 0} + userData: + assetBundleName: + assetBundleVariant: diff --git a/Unity/Plugins/SVMModel.dll b/Unity/Plugins/SVMModel.dll new file mode 100644 index 0000000..b78b629 Binary files /dev/null and b/Unity/Plugins/SVMModel.dll differ diff --git a/Unity/Plugins/SVMModel.dll.meta b/Unity/Plugins/SVMModel.dll.meta new file mode 100644 index 0000000..344264a --- /dev/null +++ b/Unity/Plugins/SVMModel.dll.meta @@ -0,0 +1,92 @@ +fileFormatVersion: 2 +guid: 930906668f1b67a42b2bb844bd098a26 +PluginImporter: + externalObjects: {} + serializedVersion: 2 + iconMap: {} + executionOrder: {} + isPreloaded: 0 + isOverridable: 0 + platformData: + - first: + '': Any + second: + enabled: 0 + settings: + Exclude Editor: 0 + Exclude Linux: 0 + Exclude Linux64: 0 + Exclude LinuxUniversal: 0 + Exclude OSXUniversal: 0 + Exclude WebGL: 0 + Exclude Win: 1 + Exclude Win64: 0 + - first: + Any: + second: + enabled: 1 + settings: {} + - first: + Editor: Editor + second: + enabled: 1 + settings: + CPU: AnyCPU + DefaultValueInitialized: true + OS: AnyOS + - first: + Facebook: Win + second: + enabled: 0 + settings: + CPU: None + - first: + Facebook: Win64 + second: + enabled: 0 + settings: + CPU: AnyCPU + - first: + Standalone: Linux + second: + enabled: 1 + settings: + CPU: x86 + - first: + Standalone: Linux64 + second: + enabled: 1 + settings: + CPU: x86_64 + - first: + Standalone: LinuxUniversal + second: + enabled: 1 + settings: + CPU: AnyCPU + - first: + Standalone: OSXUniversal + second: + enabled: 1 + settings: + CPU: AnyCPU + - first: + Standalone: Win + second: + enabled: 0 + settings: + CPU: None + - first: + Standalone: Win64 + second: + enabled: 1 + settings: + CPU: AnyCPU + - first: + WebGL: WebGL + second: + enabled: 1 + settings: {} + userData: + assetBundleName: + assetBundleVariant: