1 //
2 //  Stable_Set.cpp
3 //  Gravity
4 //
5 //  Created by Hassan Hijazi on Nov 21 2018
6 //
7 //
8 
9 
10 #include <stdio.h>
11 #include <iostream>
12 #include <string>
13 #include <stdio.h>
14 #include <cstring>
15 #include <fstream>
16 #include <gravity/Net.h>
17 #include <gravity/model.h>
18 #include <gravity/solver.h>
19 #include <DataSet.h>
20 #include <stdio.h>
21 #include <stdlib.h>
22 
23 using namespace std;
24 using namespace gravity;
25 
26 /* Builds the static SVM model */
build_svm(const DataSet<> & training_set,double mu)27 unique_ptr<Model<>> build_svm(const DataSet<>& training_set, double mu){
28 
29     /* Defining parameters ans indices */
30     auto nf = training_set._nb_features;
31     auto m1 = training_set._class_sizes[0];
32     auto m2 = training_set._class_sizes[1];
33     auto F = training_set.get_features_matrices(); /* Features matrices, one per-class */
34 
35     /* Model */
36     auto SVM = unique_ptr<Model<>>(new Model<>());
37 
38     /* Variables */
39     var<> w("w");
40     var<> xi1("xi1", pos_);
41     var<> xi2("xi2", pos_);
42     var<> b("b");
43     SVM->add(w.in(R(nf)), xi1.in(R(m1)), xi2.in(R(m2)), b.in(R(1)));
44     /* Objective function */
45     SVM->min(product(0.5,pow(w,2)) + mu*sum(xi1) + mu*sum(xi2));
46 
47     /* Constraints */
48     /* Class 1 constraints */
49     auto Class1 = Constraint<>("Class1");
50     Class1 = product(F[0],w) + b + (xi1 - 1);
51     SVM->add(Class1 >= 0);
52     /* Class 2 constraints */
53     auto Class2 = Constraint<>("Class2");
54     Class2 = product(F[1],w) + b + (1 - xi2);
55     SVM->add(Class2 <= 0);
56 
57     return SVM;
58 }
59 
60 
61 /* Builds the dual SVM model */
build_svm_dual(const DataSet<> & training_set,double mu,const string & kernel_type,double gamma,double r,unsigned d)62 unique_ptr<Model<>> build_svm_dual(const DataSet<>& training_set, double mu, const string& kernel_type, double gamma, double r, unsigned d){
63 
64     /* Defining parameters ans indices */
65     auto m = training_set._nb_points;
66     auto K = training_set.get_kernel_matrix(kernel_type, gamma, r, d);
67     auto y = training_set.get_classes();
68 
69     /* Model */
70     auto SVM = unique_ptr<Model<>>(new Model<>());
71 
72     /* Variables */
73     var<> alpha("��", 0, mu);
74     SVM->add(alpha.in(R(m)));
75 //    SVM->initialize_uniform();
76     /* Objective function */
77     SVM->min(0.5*alpha.tr()*K*alpha - sum(alpha));
78 
79     /* Constraints */
80     /* Equality constraints */
81     auto Equ0 = Constraint<>("Equation");
82     Equ0 = y.tr()*alpha;
83     SVM->add(Equ0 == 0);
84 
85     return SVM;
86 }
87 
88 
build_lazy_svm(const DataSet<> & training_set,int nb_c,double mu)89 unique_ptr<Model<>> build_lazy_svm(const DataSet<>& training_set, int nb_c, double mu){
90 
91     /* Defining parameters ans indices */
92     auto nf = training_set._nb_features;
93     auto m = training_set._nb_points;
94     auto m1 = training_set._class_sizes[0]; /* Number of points in first class */
95     auto m2 = training_set._class_sizes[1]; /* Number of points in second class */
96     auto f = training_set.get_features(); /* Feature values */
97 
98     /* Model */
99     unique_ptr<Model<>> SVM = unique_ptr<Model<>>(new Model<>());
100     /* Variables */
101     var<> w("w");
102     var<> xi("xi", pos_);
103     var<> b("b");
104     SVM->add(w.in(R(nf)), xi.in(R(m)), b.in(R(1)));
105     /* Objective function */
106     SVM->min(product(0.5,pow(w,2)) + mu*sum(xi));
107 
108     /* Constraints */
109     size_t nb_c1 = m1;
110     if (nb_c >= 0 && nb_c <= m1) { /* Activating constraint generation */
111         nb_c1 = nb_c;
112     }
113     for (auto i = 0; i<nb_c1; i++) {
114         auto Class1 = Constraint<>("Class1_"+to_string(i));
115         Class1 = product(f[0][i],w) + b + (xi(i) - 1);
116         SVM->add(Class1 >= 0);
117     }
118     /* Remaining constraints are lazy */
119     for (auto i = nb_c1; i<m1; i++) {
120         auto Class1 = Constraint<>("Class1_"+to_string(i));
121         Class1 = product(f[0][i],w) + b + (xi(i) - 1);
122         SVM->add_lazy(Class1 >= 0);
123     }
124 
125     /* Class 2 constraints */
126     size_t nb_c2 = m2;
127     if (nb_c >= 0 && nb_c < m2) { /* Activating constraint generation */
128         nb_c2 = nb_c;
129     }
130     for (auto i = 0; i<nb_c2; i++) {
131         auto Class2 = Constraint<>("Class2_"+to_string(i));
132         Class2 = product(f[1][i],w) + b + (1 - xi(i+m1));
133         SVM->add(Class2 <= 0);
134     }
135     /* Remaining constraints are lazy */
136     for (auto i = nb_c2; i<m2; i++) {
137         auto Class2 = Constraint<>("Class2_"+to_string(i));
138         Class2 = product(f[1][i],w) + b + (1 - xi(i+m1));
139         SVM->add_lazy(Class2 <= 0);
140     }
141     return SVM;
142 }
143 
144 
main(int argc,char * argv[])145 int main (int argc, char * argv[])
146 {
147     double tol = 1e-3, r = 0, gamma = 0, mu;
148     unsigned d = 3;
149     int nb_c, output = 5;
150     bool dual = false;
151     auto total_time_start = get_wall_time();
152     string solver_str ="ipopt", dual_str = "no", lazy_str = "no", mu_str = "1e+6", nb_c_str = "200", output_str = "5", kernel = "linear";
153     std::cout << "WELCOME, THIS IS AN IMPLEMENTATION OF A SUPPORT VECTOR MACHINE IN GRAVITY\n";
154 
155     string fname = string(prj_dir)+"/data_sets/Classification/Archive/svmguide1";
156     /* Reading input data */
157     DataSet<> training_set;
158     training_set.parse(fname);
159     training_set.print_stats();
160     auto nf = training_set._nb_features;
161     gamma = 1./nf;
162     SolverType solv_type = ipopt;
163 #ifdef USE_OPT_PARSER
164     /* Create a OptionParser with options */
165     op::OptionParser opt;
166     opt.add_option("h", "help",
167                    "shows option help"); // no default value means boolean options, which default value is false
168     opt.add_option("f", "file", "Input file name", fname);
169     opt.add_option("s", "solver", "Solvers: Ipopt/Cplex/Gurobi, default = Ipopt", solver_str);
170     opt.add_option("d", "dual", "Solve dual SVM, yes/no, default = no", dual_str);
171     opt.add_option("k", "kernel", "Choose kernel, linear/poly/rbf/sigm, default = linear", kernel);
172     opt.add_option("o", "output", "Output level, default = 5", output_str);
173     opt.add_option("lz", "lazy", "Generate constraints in a lazy fashion, yes/no, default = no", lazy_str);
174     opt.add_option("mu", "multiplier", "Value of penalization multiplier, default = 1e+6", mu_str);
175     opt.add_option("nb", "nb_init_cstr", "Initial number of constraints to add, default = 200 (will add the constraints corresponding to the first 'nb' points in each class). Setting this parameter to -1 will deactivate constraint generation.", nb_c_str);
176     /* Parse the options and verify that all went well. If not, errors and help will be shown */
177     bool correct_parsing = opt.parse_options(argc, argv);
178 
179     if (!correct_parsing) {
180         return EXIT_FAILURE;
181     }
182 
183     fname = opt["f"];
184     bool has_help = op::str2bool(opt["h"]);
185     if (has_help) {
186         opt.show_help();
187         exit(0);
188     }
189     solver_str = opt["s"];
190     if (solver_str.compare("Gurobi")==0) {
191         solv_type = gurobi;
192     }
193     else if(solver_str.compare("Cplex")==0) {
194         solv_type = cplex;
195     }else if(solver_str.compare("Mosek")==0) {
196         solv_type = _mosek;
197     }
198     lazy_str = opt["lz"];
199     bool lazy = false;
200     if (lazy_str.compare("yes")==0) {
201         lazy = true;
202     }
203     dual_str = opt["d"];
204     if (dual_str.compare("yes")==0) {
205         dual = true;
206     }
207     mu_str = opt["mu"];
208     mu = op::str2double(mu_str);
209 
210     nb_c_str = opt["nb"];
211     nb_c = op::str2int(nb_c_str);
212 
213     output_str = opt["o"];
214     output = op::str2int(output_str);
215 
216     unique_ptr<Model<>> SVM;
217     if(dual){
218         kernel = opt["k"];
219         SVM = build_svm_dual(training_set, mu, kernel,1./training_set._nb_features,0,3);
220     }
221     else {
222         if (lazy && nb_c>0) {
223             SVM = build_lazy_svm(training_set, nb_c, mu);
224         }
225         else {
226             SVM = build_svm(training_set, mu);
227         }
228     }
229 #else
230     auto SVM = build_svm(training_set, mu);
231 #endif
232     SVM->print_symbolic();
233 //    SVM->print();
234     /* Start Timers */
235     solver<> SVM_solver(*SVM,solv_type);
236     double solver_time_start = get_wall_time();
237     SVM_solver.run(output,tol);
238     double solver_time_end = get_wall_time();
239     double total_time_end = get_wall_time();
240     auto solve_time = solver_time_end - solver_time_start;
241     auto total_time = total_time_end - total_time_start;
242     /* Uncomment line below to print expanded model */
243     /* SVM->print(); */
244 //    SVM->print_solution(false);
245     /* Testing */
246     unsigned err = 0;
247     DataSet<> test_set;
248     test_set.parse(fname+".t", &training_set);
249     if(dual){
250         double b = 0;
251         auto alpha = SVM->get_var<double>("��");
252         auto y = training_set.get_classes();
253         bool point_on_margin = false;
254         for (auto i = 0; i<training_set._nb_points; i++) {
255             if (alpha.eval(i)>tol && alpha.eval(i)<mu-tol) {//point lies on the margin
256                 auto pi = training_set._all_points[i];
257                 b = y.eval(i);
258                 for (auto k = 0; k<training_set._nb_points; k++) {
259                     auto pk = training_set._all_points[k];
260                     b -= alpha.eval(k)*y.eval(k)*training_set.K(*pi,*pk, kernel, gamma, r, d);
261                 }
262                 point_on_margin = true;
263                 break;
264             }
265         }
266         if (!point_on_margin) {
267             DebugOn("No points found on margin, aborting" << endl);
268             return 0;
269         }
270         DebugOn("b = " << b << endl);
271         for (auto c = 0; c<test_set._nb_classes; c++) {
272             for (auto i = 0; i<test_set._class_sizes[c]; i++) {
273                 double lhs = 0;
274                 auto pi = &test_set._points[c][i];
275                 for (auto k = 0; k<training_set._nb_points; k++) {
276                     auto pk = training_set._all_points[k];
277                     lhs += alpha.eval(k)*y.eval(k)*training_set.K(*pi,*pk, kernel, gamma, r, d);
278                 }
279                 lhs += b;
280                 if (c == 0 && lhs <= 0) {
281                     err++;
282                 }
283                 else if (c == 1 && lhs > 0) {
284                     err++;
285                 }
286             }
287         }
288     }
289     else {
290         auto w = SVM->get_var<double>("w");
291         auto b = SVM->get_var<double>("b");
292         DebugOn("b = " << b.eval() << endl);
293         for (auto c = 0; c<test_set._nb_classes; c++) {
294             for (auto i = 0; i<test_set._class_sizes[c]; i++) {
295                 double lhs = 0;
296                 auto pt = test_set._points[c][i];
297                 for (auto j = 0; j<nf; j++) {
298                     lhs += pt._features[j]*w.eval(j);
299                 }
300                 lhs += b.eval();
301                 if (c == 0 && lhs <= 0) {
302                     err++;
303                 }
304                 else if (c == 1 && lhs > 0) {
305                     err++;
306                 }
307             }
308         }
309 
310     }
311     cout << "Number of misclassified = " << err << " out of " << test_set._nb_points << endl;
312     cout << "Accuracy = " << 100.*(test_set._nb_points-err)/test_set._nb_points << "%" << endl;
313 
314     /** Terminal output */
315     cout << "Done running the SVM model\n";
316     cout << "Solve wall clock time =  " << solve_time << "\n";
317     cout << "Total wall clock time =  " << total_time << "\n";
318     return 0;
319 }
320 
321 
322 
323 
324