Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.LibSVM/1.6.3/LibSVM-1.6.3/Solver.cs @ 3884

Last change on this file since 3884 was 3884, checked in by gkronber, 14 years ago

Worked on support vector regression operators and views. #1009

File size: 79.2 KB
Line 
1/*
2 * SVM.NET Library
3 * Copyright (C) 2008 Matthew Johnson
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
17 */
18
19
20using System;
21using System.Linq;
22using System.Collections.Generic;
23using System.Diagnostics;
24using System.IO;
25
26namespace SVM
27{
28
29    // An SMO algorithm in Fan et al., JMLR 6(2005), p. 1889--1918
30    // Solves:
31    //
32    //  Min 0.5(\alpha^T Q \alpha) + p^T \alpha
33    //
34    //    y^T \alpha = \delta
35    //    y_i = +1 or -1
36    //    0 <= alpha_i <= Cp for y_i = 1
37    //    0 <= alpha_i <= Cn for y_i = -1
38    //
39    // Given:
40    //
41    //  Q, p, y, Cp, Cn, and an initial feasible point \alpha
42    //  l is the size of vectors and matrices
43    //  eps is the stopping tolerance
44    //
45    // solution will be put in \alpha, objective value will be put in obj
46    //
47    internal class Solver
48    {
49        protected int active_size;
50        protected sbyte[] y;
51        protected double[] G;   // gradient of objective function
52        private const byte LOWER_BOUND = 0;
53        private const byte UPPER_BOUND = 1;
54        private const byte FREE = 2;
55        private byte[] alpha_status;  // LOWER_BOUND, UPPER_BOUND, FREE
56        private double[] alpha;
57        protected IQMatrix Q;
58        protected float[] QD;
59        protected double EPS;
60        private double Cp, Cn;
61        private double[] p;
62        private int[] active_set;
63        private double[] G_bar;   // gradient, if we treat free variables as 0
64        protected int l;
65        protected bool unshrink;  // XXX
66
67        protected const double INF = double.PositiveInfinity;
68
69        private double get_C(int i)
70        {
71            return (y[i] > 0) ? Cp : Cn;
72        }
73
74        private void update_alpha_status(int i)
75        {
76            if (alpha[i] >= get_C(i))
77                alpha_status[i] = UPPER_BOUND;
78            else if (alpha[i] <= 0)
79                alpha_status[i] = LOWER_BOUND;
80            else alpha_status[i] = FREE;
81        }
82
83        protected bool is_upper_bound(int i) { return alpha_status[i] == UPPER_BOUND; }
84        protected bool is_lower_bound(int i) { return alpha_status[i] == LOWER_BOUND; }
85
86        private bool is_free(int i) { return alpha_status[i] == FREE; }
87
88        public class SolutionInfo
89        {
90            public double obj;
91            public double rho;
92            public double upper_bound_p;
93            public double upper_bound_n;
94            public double r;  // for Solver_NU
95        }
96
97        protected void swap_index(int i, int j)
98        {
99            Q.SwapIndex(i, j);
100            y.SwapIndex(i, j);
101            G.SwapIndex(i, j);
102            alpha_status.SwapIndex(i, j);
103            alpha.SwapIndex(i, j);
104            p.SwapIndex(i, j);
105            active_set.SwapIndex(i, j);
106            G_bar.SwapIndex(i, j);
107        }
108
109        protected void reconstruct_gradient()
110        {
111            // reconstruct inactive elements of G from G_bar and free variables
112
113            if (active_size == l) return;
114
115            int i, j;
116            int nr_free = 0;
117
118            for (j = active_size; j < l; j++)
119                G[j] = G_bar[j] + p[j];
120
121            for (j = 0; j < active_size; j++)
122                if (is_free(j))
123                    nr_free++;
124
125            if (2 * nr_free < active_size)
126                Procedures.info("\nWarning: using -h 0 may be faster\n");
127
128            if (nr_free * l > 2 * active_size * (l - active_size))
129            {
130                for (i = active_size; i < l; i++)
131                {
132                    float[] Q_i = Q.GetQ(i, active_size);
133                    for (j = 0; j < active_size; j++)
134                        if (is_free(j))
135                            G[i] += alpha[j] * Q_i[j];
136                }
137            }
138            else
139            {
140                for (i = 0; i < active_size; i++)
141                    if (is_free(i))
142                    {
143                        float[] Q_i = Q.GetQ(i, l);
144                        double alpha_i = alpha[i];
145                        for (j = active_size; j < l; j++)
146                            G[j] += alpha_i * Q_i[j];
147                    }
148            }
149        }
150
151        public virtual void Solve(int l, IQMatrix Q, double[] p_, sbyte[] y_, double[] alpha_, double Cp, double Cn, double eps, SolutionInfo si, bool shrinking)
152        {
153            this.l = l;
154            this.Q = Q;
155            QD = Q.GetQD();
156            p = (double[])p_.Clone();
157            y = (sbyte[])y_.Clone();
158            alpha = (double[])alpha_.Clone();
159            this.Cp = Cp;
160            this.Cn = Cn;
161            this.EPS = eps;
162            this.unshrink = false;
163
164            // initialize alpha_status
165            {
166                alpha_status = new byte[l];
167                for (int i = 0; i < l; i++)
168                    update_alpha_status(i);
169            }
170
171            // initialize active set (for shrinking)
172            {
173                active_set = new int[l];
174                for (int i = 0; i < l; i++)
175                    active_set[i] = i;
176                active_size = l;
177            }
178
179            // initialize gradient
180            {
181                G = new double[l];
182                G_bar = new double[l];
183                int i;
184                for (i = 0; i < l; i++)
185                {
186                    G[i] = p[i];
187                    G_bar[i] = 0;
188                }
189                for (i = 0; i < l; i++)
190                    if (!is_lower_bound(i))
191                    {
192                        float[] Q_i = Q.GetQ(i, l);
193                        double alpha_i = alpha[i];
194                        int j;
195                        for (j = 0; j < l; j++)
196                            G[j] += alpha_i * Q_i[j];
197                        if (is_upper_bound(i))
198                            for (j = 0; j < l; j++)
199                                G_bar[j] += get_C(i) * Q_i[j];
200                    }
201            }
202
203            // optimization step
204
205            int iter = 0;
206            int counter = Math.Min(l, 1000) + 1;
207            int[] working_set = new int[2];
208
209            while (true)
210            {
211                // show progress and do shrinking
212
213                if (--counter == 0)
214                {
215                    counter = Math.Min(l, 1000);
216                    if (shrinking) do_shrinking();
217                    Procedures.info(".");
218                }
219
220                if (select_working_set(working_set) != 0)
221                {
222                    // reconstruct the whole gradient
223                    reconstruct_gradient();
224                    // reset active set size and check
225                    active_size = l;
226                    Procedures.info("*");
227                    if (select_working_set(working_set) != 0)
228                        break;
229                    else
230                        counter = 1;  // do shrinking next iteration
231                }
232
233                int i = working_set[0];
234                int j = working_set[1];
235
236                ++iter;
237
238                // update alpha[i] and alpha[j], handle bounds carefully
239
240                float[] Q_i = Q.GetQ(i, active_size);
241                float[] Q_j = Q.GetQ(j, active_size);
242
243                double C_i = get_C(i);
244                double C_j = get_C(j);
245
246                double old_alpha_i = alpha[i];
247                double old_alpha_j = alpha[j];
248
249                if (y[i] != y[j])
250                {
251                    double quad_coef = Q_i[i] + Q_j[j] + 2 * Q_i[j];
252                    if (quad_coef <= 0)
253                        quad_coef = 1e-12;
254                    double delta = (-G[i] - G[j]) / quad_coef;
255                    double diff = alpha[i] - alpha[j];
256                    alpha[i] += delta;
257                    alpha[j] += delta;
258
259                    if (diff > 0)
260                    {
261                        if (alpha[j] < 0)
262                        {
263                            alpha[j] = 0;
264                            alpha[i] = diff;
265                        }
266                    }
267                    else
268                    {
269                        if (alpha[i] < 0)
270                        {
271                            alpha[i] = 0;
272                            alpha[j] = -diff;
273                        }
274                    }
275                    if (diff > C_i - C_j)
276                    {
277                        if (alpha[i] > C_i)
278                        {
279                            alpha[i] = C_i;
280                            alpha[j] = C_i - diff;
281                        }
282                    }
283                    else
284                    {
285                        if (alpha[j] > C_j)
286                        {
287                            alpha[j] = C_j;
288                            alpha[i] = C_j + diff;
289                        }
290                    }
291                }
292                else
293                {
294                    double quad_coef = Q_i[i] + Q_j[j] - 2 * Q_i[j];
295                    if (quad_coef <= 0)
296                        quad_coef = 1e-12;
297                    double delta = (G[i] - G[j]) / quad_coef;
298                    double sum = alpha[i] + alpha[j];
299                    alpha[i] -= delta;
300                    alpha[j] += delta;
301
302                    if (sum > C_i)
303                    {
304                        if (alpha[i] > C_i)
305                        {
306                            alpha[i] = C_i;
307                            alpha[j] = sum - C_i;
308                        }
309                    }
310                    else
311                    {
312                        if (alpha[j] < 0)
313                        {
314                            alpha[j] = 0;
315                            alpha[i] = sum;
316                        }
317                    }
318                    if (sum > C_j)
319                    {
320                        if (alpha[j] > C_j)
321                        {
322                            alpha[j] = C_j;
323                            alpha[i] = sum - C_j;
324                        }
325                    }
326                    else
327                    {
328                        if (alpha[i] < 0)
329                        {
330                            alpha[i] = 0;
331                            alpha[j] = sum;
332                        }
333                    }
334                }
335
336                // update G
337
338                double delta_alpha_i = alpha[i] - old_alpha_i;
339                double delta_alpha_j = alpha[j] - old_alpha_j;
340
341                for (int k = 0; k < active_size; k++)
342                {
343                    G[k] += Q_i[k] * delta_alpha_i + Q_j[k] * delta_alpha_j;
344                }
345
346                // update alpha_status and G_bar
347
348                {
349                    bool ui = is_upper_bound(i);
350                    bool uj = is_upper_bound(j);
351                    update_alpha_status(i);
352                    update_alpha_status(j);
353                    int k;
354                    if (ui != is_upper_bound(i))
355                    {
356                        Q_i = Q.GetQ(i, l);
357                        if (ui)
358                            for (k = 0; k < l; k++)
359                                G_bar[k] -= C_i * Q_i[k];
360                        else
361                            for (k = 0; k < l; k++)
362                                G_bar[k] += C_i * Q_i[k];
363                    }
364
365                    if (uj != is_upper_bound(j))
366                    {
367                        Q_j = Q.GetQ(j, l);
368                        if (uj)
369                            for (k = 0; k < l; k++)
370                                G_bar[k] -= C_j * Q_j[k];
371                        else
372                            for (k = 0; k < l; k++)
373                                G_bar[k] += C_j * Q_j[k];
374                    }
375                }
376
377            }
378
379            // calculate rho
380
381            si.rho = calculate_rho();
382
383            // calculate objective value
384            {
385                double v = 0;
386                int i;
387                for (i = 0; i < l; i++)
388                    v += alpha[i] * (G[i] + p[i]);
389
390                si.obj = v / 2;
391            }
392
393            // put back the solution
394            {
395                for (int i = 0; i < l; i++)
396                    alpha_[active_set[i]] = alpha[i];
397            }
398
399            si.upper_bound_p = Cp;
400            si.upper_bound_n = Cn;
401
402            Procedures.info("\noptimization finished, #iter = " + iter + "\n");
403        }
404
405        // return 1 if already optimal, return 0 otherwise
406        protected virtual int select_working_set(int[] working_set)
407        {
408            // return i,j such that
409            // i: Maximizes -y_i * grad(f)_i, i in I_up(\alpha)
410            // j: mimimizes the decrease of obj value
411            //    (if quadratic coefficeint <= 0, replace it with tau)
412            //    -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha)
413
414            double GMax = -INF;
415            double GMax2 = -INF;
416            int GMax_idx = -1;
417            int GMin_idx = -1;
418            double obj_diff_Min = INF;
419
420            for (int t = 0; t < active_size; t++)
421                if (y[t] == +1)
422                {
423                    if (!is_upper_bound(t))
424                        if (-G[t] >= GMax)
425                        {
426                            GMax = -G[t];
427                            GMax_idx = t;
428                        }
429                }
430                else
431                {
432                    if (!is_lower_bound(t))
433                        if (G[t] >= GMax)
434                        {
435                            GMax = G[t];
436                            GMax_idx = t;
437                        }
438                }
439
440            int i = GMax_idx;
441            float[] Q_i = null;
442            if (i != -1) // null Q_i not accessed: GMax=-INF if i=-1
443                Q_i = Q.GetQ(i, active_size);
444
445            for (int j = 0; j < active_size; j++)
446            {
447                if (y[j] == +1)
448                {
449                    if (!is_lower_bound(j))
450                    {
451                        double grad_diff = GMax + G[j];
452                        if (G[j] >= GMax2)
453                            GMax2 = G[j];
454                        if (grad_diff > 0)
455                        {
456                            double obj_diff;
457                            double quad_coef = Q_i[i] + QD[j] - 2.0 * y[i] * Q_i[j];
458                            if (quad_coef > 0)
459                                obj_diff = -(grad_diff * grad_diff) / quad_coef;
460                            else
461                                obj_diff = -(grad_diff * grad_diff) / 1e-12;
462
463                            if (obj_diff <= obj_diff_Min)
464                            {
465                                GMin_idx = j;
466                                obj_diff_Min = obj_diff;
467                            }
468                        }
469                    }
470                }
471                else
472                {
473                    if (!is_upper_bound(j))
474                    {
475                        double grad_diff = GMax - G[j];
476                        if (-G[j] >= GMax2)
477                            GMax2 = -G[j];
478                        if (grad_diff > 0)
479                        {
480                            double obj_diff;
481                            double quad_coef = Q_i[i] + QD[j] + 2.0 * y[i] * Q_i[j];
482                            if (quad_coef > 0)
483                                obj_diff = -(grad_diff * grad_diff) / quad_coef;
484                            else
485                                obj_diff = -(grad_diff * grad_diff) / 1e-12;
486
487                            if (obj_diff <= obj_diff_Min)
488                            {
489                                GMin_idx = j;
490                                obj_diff_Min = obj_diff;
491                            }
492                        }
493                    }
494                }
495            }
496
497            if (GMax + GMax2 < EPS)
498                return 1;
499
500            working_set[0] = GMax_idx;
501            working_set[1] = GMin_idx;
502            return 0;
503        }
504
505        private bool be_shrunk(int i, double GMax1, double GMax2)
506        {
507            if (is_upper_bound(i))
508            {
509                if (y[i] == +1)
510                    return (-G[i] > GMax1);
511                else
512                    return (-G[i] > GMax2);
513            }
514            else if (is_lower_bound(i))
515            {
516                if (y[i] == +1)
517                    return (G[i] > GMax2);
518                else
519                    return (G[i] > GMax1);
520            }
521            else
522                return (false);
523        }
524
525        protected virtual void do_shrinking()
526        {
527            int i;
528            double GMax1 = -INF;    // Max { -y_i * grad(f)_i | i in I_up(\alpha) }
529            double GMax2 = -INF;    // Max { y_i * grad(f)_i | i in I_low(\alpha) }
530
531            // find Maximal violating pair first
532            for (i = 0; i < active_size; i++)
533            {
534                if (y[i] == +1)
535                {
536                    if (!is_upper_bound(i))
537                    {
538                        if (-G[i] >= GMax1)
539                            GMax1 = -G[i];
540                    }
541                    if (!is_lower_bound(i))
542                    {
543                        if (G[i] >= GMax2)
544                            GMax2 = G[i];
545                    }
546                }
547                else
548                {
549                    if (!is_upper_bound(i))
550                    {
551                        if (-G[i] >= GMax2)
552                            GMax2 = -G[i];
553                    }
554                    if (!is_lower_bound(i))
555                    {
556                        if (G[i] >= GMax1)
557                            GMax1 = G[i];
558                    }
559                }
560            }
561
562            if (unshrink == false && GMax1 + GMax2 <= EPS * 10)
563            {
564                unshrink = true;
565                reconstruct_gradient();
566                active_size = l;
567            }
568
569            for (i = 0; i < active_size; i++)
570                if (be_shrunk(i, GMax1, GMax2))
571                {
572                    active_size--;
573                    while (active_size > i)
574                    {
575                        if (!be_shrunk(active_size, GMax1, GMax2))
576                        {
577                            swap_index(i, active_size);
578                            break;
579                        }
580                        active_size--;
581                    }
582                }
583        }
584
585        protected virtual double calculate_rho()
586        {
587            double r;
588            int nr_free = 0;
589            double ub = INF, lb = -INF, sum_free = 0;
590            for (int i = 0; i < active_size; i++)
591            {
592                double yG = y[i] * G[i];
593
594                if (is_lower_bound(i))
595                {
596                    if (y[i] > 0)
597                        ub = Math.Min(ub, yG);
598                    else
599                        lb = Math.Max(lb, yG);
600                }
601                else if (is_upper_bound(i))
602                {
603                    if (y[i] < 0)
604                        ub = Math.Min(ub, yG);
605                    else
606                        lb = Math.Max(lb, yG);
607                }
608                else
609                {
610                    ++nr_free;
611                    sum_free += yG;
612                }
613            }
614
615            if (nr_free > 0)
616                r = sum_free / nr_free;
617            else
618                r = (ub + lb) / 2;
619
620            return r;
621        }
622
623    }
624
625    //
626    // Solver for nu-svm classification and regression
627    //
628    // additional constraint: e^T \alpha = constant
629    //
630    class Solver_NU : Solver
631    {
632        private SolutionInfo si;
633
634        public sealed override void Solve(int l, IQMatrix Q, double[] p, sbyte[] y,
635               double[] alpha, double Cp, double Cn, double eps,
636               SolutionInfo si, bool shrinking)
637        {
638            this.si = si;
639            base.Solve(l, Q, p, y, alpha, Cp, Cn, eps, si, shrinking);
640        }
641
642        // return 1 if already optimal, return 0 otherwise
643        protected override sealed int select_working_set(int[] working_set)
644        {
645            // return i,j such that y_i = y_j and
646            // i: Maximizes -y_i * grad(f)_i, i in I_up(\alpha)
647            // j: Minimizes the decrease of obj value
648            //    (if quadratic coefficeint <= 0, replace it with tau)
649            //    -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha)
650
651            double GMaxp = -INF;
652            double GMaxp2 = -INF;
653            int GMaxp_idx = -1;
654
655            double GMaxn = -INF;
656            double GMaxn2 = -INF;
657            int GMaxn_idx = -1;
658
659            int GMin_idx = -1;
660            double obj_diff_Min = INF;
661
662            for (int t = 0; t < active_size; t++)
663                if (y[t] == +1)
664                {
665                    if (!is_upper_bound(t))
666                        if (-G[t] >= GMaxp)
667                        {
668                            GMaxp = -G[t];
669                            GMaxp_idx = t;
670                        }
671                }
672                else
673                {
674                    if (!is_lower_bound(t))
675                        if (G[t] >= GMaxn)
676                        {
677                            GMaxn = G[t];
678                            GMaxn_idx = t;
679                        }
680                }
681
682            int ip = GMaxp_idx;
683            int iN = GMaxn_idx;
684            float[] Q_ip = null;
685            float[] Q_in = null;
686            if (ip != -1) // null Q_ip not accessed: GMaxp=-INF if ip=-1
687                Q_ip = Q.GetQ(ip, active_size);
688            if (iN != -1)
689                Q_in = Q.GetQ(iN, active_size);
690
691            for (int j = 0; j < active_size; j++)
692            {
693                if (y[j] == +1)
694                {
695                    if (!is_lower_bound(j))
696                    {
697                        double grad_diff = GMaxp + G[j];
698                        if (G[j] >= GMaxp2)
699                            GMaxp2 = G[j];
700                        if (grad_diff > 0)
701                        {
702                            double obj_diff;
703                            double quad_coef = Q_ip[ip] + QD[j] - 2 * Q_ip[j];
704                            if (quad_coef > 0)
705                                obj_diff = -(grad_diff * grad_diff) / quad_coef;
706                            else
707                                obj_diff = -(grad_diff * grad_diff) / 1e-12;
708
709                            if (obj_diff <= obj_diff_Min)
710                            {
711                                GMin_idx = j;
712                                obj_diff_Min = obj_diff;
713                            }
714                        }
715                    }
716                }
717                else
718                {
719                    if (!is_upper_bound(j))
720                    {
721                        double grad_diff = GMaxn - G[j];
722                        if (-G[j] >= GMaxn2)
723                            GMaxn2 = -G[j];
724                        if (grad_diff > 0)
725                        {
726                            double obj_diff;
727                            double quad_coef = Q_in[iN] + QD[j] - 2 * Q_in[j];
728                            if (quad_coef > 0)
729                                obj_diff = -(grad_diff * grad_diff) / quad_coef;
730                            else
731                                obj_diff = -(grad_diff * grad_diff) / 1e-12;
732
733                            if (obj_diff <= obj_diff_Min)
734                            {
735                                GMin_idx = j;
736                                obj_diff_Min = obj_diff;
737                            }
738                        }
739                    }
740                }
741            }
742
743            if (Math.Max(GMaxp + GMaxp2, GMaxn + GMaxn2) < EPS)
744                return 1;
745
746            if (y[GMin_idx] == +1)
747                working_set[0] = GMaxp_idx;
748            else
749                working_set[0] = GMaxn_idx;
750            working_set[1] = GMin_idx;
751
752            return 0;
753        }
754
755        private bool be_shrunk(int i, double GMax1, double GMax2, double GMax3, double GMax4)
756        {
757            if (is_upper_bound(i))
758            {
759                if (y[i] == +1)
760                    return (-G[i] > GMax1);
761                else
762                    return (-G[i] > GMax4);
763            }
764            else if (is_lower_bound(i))
765            {
766                if (y[i] == +1)
767                    return (G[i] > GMax2);
768                else
769                    return (G[i] > GMax3);
770            }
771            else
772                return (false);
773        }
774
775        protected override sealed void do_shrinking()
776        {
777            double GMax1 = -INF;  // Max { -y_i * grad(f)_i | y_i = +1, i in I_up(\alpha) }
778            double GMax2 = -INF;  // Max { y_i * grad(f)_i | y_i = +1, i in I_low(\alpha) }
779            double GMax3 = -INF;  // Max { -y_i * grad(f)_i | y_i = -1, i in I_up(\alpha) }
780            double GMax4 = -INF;  // Max { y_i * grad(f)_i | y_i = -1, i in I_low(\alpha) }
781
782            // find Maximal violating pair first
783            int i;
784            for (i = 0; i < active_size; i++)
785            {
786                if (!is_upper_bound(i))
787                {
788                    if (y[i] == +1)
789                    {
790                        if (-G[i] > GMax1) GMax1 = -G[i];
791                    }
792                    else if (-G[i] > GMax4) GMax4 = -G[i];
793                }
794                if (!is_lower_bound(i))
795                {
796                    if (y[i] == +1)
797                    {
798                        if (G[i] > GMax2) GMax2 = G[i];
799                    }
800                    else if (G[i] > GMax3) GMax3 = G[i];
801                }
802            }
803
804            if (unshrink == false && Math.Max(GMax1 + GMax2, GMax3 + GMax4) <= EPS * 10)
805            {
806                unshrink = true;
807                reconstruct_gradient();
808                active_size = l;
809            }
810
811            for (i = 0; i < active_size; i++)
812                if (be_shrunk(i, GMax1, GMax2, GMax3, GMax4))
813                {
814                    active_size--;
815                    while (active_size > i)
816                    {
817                        if (!be_shrunk(active_size, GMax1, GMax2, GMax3, GMax4))
818                        {
819                            swap_index(i, active_size);
820                            break;
821                        }
822                        active_size--;
823                    }
824                }
825        }
826
827        protected override sealed double calculate_rho()
828        {
829            int nr_free1 = 0, nr_free2 = 0;
830            double ub1 = INF, ub2 = INF;
831            double lb1 = -INF, lb2 = -INF;
832            double sum_free1 = 0, sum_free2 = 0;
833
834            for (int i = 0; i < active_size; i++)
835            {
836                if (y[i] == +1)
837                {
838                    if (is_lower_bound(i))
839                        ub1 = Math.Min(ub1, G[i]);
840                    else if (is_upper_bound(i))
841                        lb1 = Math.Max(lb1, G[i]);
842                    else
843                    {
844                        ++nr_free1;
845                        sum_free1 += G[i];
846                    }
847                }
848                else
849                {
850                    if (is_lower_bound(i))
851                        ub2 = Math.Min(ub2, G[i]);
852                    else if (is_upper_bound(i))
853                        lb2 = Math.Max(lb2, G[i]);
854                    else
855                    {
856                        ++nr_free2;
857                        sum_free2 += G[i];
858                    }
859                }
860            }
861
862            double r1, r2;
863            if (nr_free1 > 0)
864                r1 = sum_free1 / nr_free1;
865            else
866                r1 = (ub1 + lb1) / 2;
867
868            if (nr_free2 > 0)
869                r2 = sum_free2 / nr_free2;
870            else
871                r2 = (ub2 + lb2) / 2;
872
873            si.r = (r1 + r2) / 2;
874            return (r1 - r2) / 2;
875        }
876    }
877
878    //
879    // Q matrices for various formulations
880    //
881    class SVC_Q : Kernel
882    {
883        private sbyte[] y;
884        private Cache cache;
885        private float[] QD;
886
887        public SVC_Q(Problem prob, Parameter param, sbyte[] y_) : base(prob.Count, prob.X, param)
888        {
889            y = (sbyte[])y_.Clone();
890            cache = new Cache(prob.Count, (long)(param.CacheSize * (1 << 20)));
891            QD = new float[prob.Count];
892            for (int i = 0; i < prob.Count; i++)
893                QD[i] = (float)KernelFunction(i, i);
894        }
895
896        public override sealed float[] GetQ(int i, int len)
897        {
898            float[] data = null;
899            int start, j;
900            if ((start = cache.GetData(i, ref data, len)) < len)
901            {
902                for (j = start; j < len; j++)
903                    data[j] = (float)(y[i] * y[j] * KernelFunction(i, j));
904            }
905            return data;
906        }
907
908        public override sealed float[] GetQD()
909        {
910            return QD;
911        }
912
913        public override sealed void SwapIndex(int i, int j)
914        {
915            cache.SwapIndex(i, j);
916            base.SwapIndex(i, j);
917            y.SwapIndex(i, j);
918            QD.SwapIndex(i, j);
919        }
920    }
921
922    class ONE_CLASS_Q : Kernel
923    {
924        private Cache cache;
925        private float[] QD;
926
927        public ONE_CLASS_Q(Problem prob, Parameter param) :  base(prob.Count, prob.X, param)
928        {
929            cache = new Cache(prob.Count, (long)(param.CacheSize * (1 << 20)));
930            QD = new float[prob.Count];
931            for (int i = 0; i < prob.Count; i++)
932                QD[i] = (float)KernelFunction(i, i);
933        }
934
935        public override sealed float[] GetQ(int i, int len)
936        {
937            float[] data = null;
938            int start, j;
939            if ((start = cache.GetData(i, ref data, len)) < len)
940            {
941                for (j = start; j < len; j++)
942                    data[j] = (float)KernelFunction(i, j);
943            }
944            return data;
945        }
946
947        public override sealed float[] GetQD()
948        {
949            return QD;
950        }
951
952        public override sealed void SwapIndex(int i, int j)
953        {
954            cache.SwapIndex(i, j);
955            base.SwapIndex(i, j);
956            QD.SwapIndex(i, j);
957        }
958    }
959
960    class SVR_Q : Kernel
961    {
962        private int l;
963        private Cache cache;
964        private sbyte[] sign;
965        private int[] index;
966        private int next_buffer;
967        private float[][] buffer;
968        private float[] QD;
969
970        public SVR_Q(Problem prob, Parameter param) : base(prob.Count, prob.X, param)
971        {
972            l = prob.Count;
973            cache = new Cache(l, (long)(param.CacheSize * (1 << 20)));
974            QD = new float[2 * l];
975            sign = new sbyte[2 * l];
976            index = new int[2 * l];
977            for (int k = 0; k < l; k++)
978            {
979                sign[k] = 1;
980                sign[k + l] = -1;
981                index[k] = k;
982                index[k + l] = k;
983                QD[k] = (float)KernelFunction(k, k);
984                QD[k + l] = QD[k];
985            }
986            buffer = new float[2][];
987            buffer[0] = new float[2 * l];
988            buffer[1] = new float[2 * l];
989            next_buffer = 0;
990        }
991
992        public override sealed void SwapIndex(int i, int j)
993        {
994            sign.SwapIndex(i, j);
995            index.SwapIndex(i, j);
996            QD.SwapIndex(i, j);
997        }
998
999        public override sealed float[] GetQ(int i, int len)
1000        {
1001            float[] data = null;
1002            int j, real_i = index[i];
1003            if (cache.GetData(real_i, ref data, l) < l)
1004            {
1005                for (j = 0; j < l; j++)
1006                    data[j] = (float)KernelFunction(real_i, j);
1007            }
1008
1009            // reorder and copy
1010            float[] buf = buffer[next_buffer];
1011            next_buffer = 1 - next_buffer;
1012            sbyte si = sign[i];
1013            for (j = 0; j < len; j++)
1014                buf[j] = (float)si * sign[j] * data[index[j]];
1015            return buf;
1016        }
1017
1018        public override sealed float[] GetQD()
1019        {
1020            return QD;
1021        }
1022    }
1023
1024    internal class Procedures
1025    {
1026        private static bool _verbose;
1027        public static bool IsVerbose
1028        {
1029            get
1030            {
1031                return _verbose;
1032            }
1033            set
1034            {
1035                _verbose = value;
1036            }
1037        }
1038        //
1039        // construct and solve various formulations
1040        //
1041        public const int LIBSVM_VERSION = 289;
1042
1043        public static TextWriter svm_print_string = Console.Out;
1044
1045        public static void info(string s)
1046        {
1047            if(_verbose)
1048                svm_print_string.Write(s);
1049        }
1050
1051        private static void solve_c_svc(Problem prob, Parameter param,
1052                        double[] alpha, Solver.SolutionInfo si,
1053                        double Cp, double Cn)
1054        {
1055            int l = prob.Count;
1056            double[] Minus_ones = new double[l];
1057            sbyte[] y = new sbyte[l];
1058
1059            int i;
1060
1061            for (i = 0; i < l; i++)
1062            {
1063                alpha[i] = 0;
1064                Minus_ones[i] = -1;
1065                if (prob.Y[i] > 0) y[i] = +1; else y[i] = -1;
1066            }
1067
1068            Solver s = new Solver();
1069            s.Solve(l, new SVC_Q(prob, param, y), Minus_ones, y,
1070                alpha, Cp, Cn, param.EPS, si, param.Shrinking);
1071
1072            double sum_alpha = 0;
1073            for (i = 0; i < l; i++)
1074                sum_alpha += alpha[i];
1075
1076            if (Cp == Cn)
1077                Procedures.info("nu = " + sum_alpha / (Cp * prob.Count) + "\n");
1078
1079            for (i = 0; i < l; i++)
1080                alpha[i] *= y[i];
1081        }
1082
1083        private static void solve_nu_svc(Problem prob, Parameter param,
1084                        double[] alpha, Solver.SolutionInfo si)
1085        {
1086            int i;
1087            int l = prob.Count;
1088            double nu = param.Nu;
1089
1090            sbyte[] y = new sbyte[l];
1091
1092            for (i = 0; i < l; i++)
1093                if (prob.Y[i] > 0)
1094                    y[i] = +1;
1095                else
1096                    y[i] = -1;
1097
1098            double sum_pos = nu * l / 2;
1099            double sum_neg = nu * l / 2;
1100
1101            for (i = 0; i < l; i++)
1102                if (y[i] == +1)
1103                {
1104                    alpha[i] = Math.Min(1.0, sum_pos);
1105                    sum_pos -= alpha[i];
1106                }
1107                else
1108                {
1109                    alpha[i] = Math.Min(1.0, sum_neg);
1110                    sum_neg -= alpha[i];
1111                }
1112
1113            double[] zeros = new double[l];
1114
1115            for (i = 0; i < l; i++)
1116                zeros[i] = 0;
1117
1118            Solver_NU s = new Solver_NU();
1119            s.Solve(l, new SVC_Q(prob, param, y), zeros, y, alpha, 1.0, 1.0, param.EPS, si, param.Shrinking);
1120            double r = si.r;
1121
1122            Procedures.info("C = " + 1 / r + "\n");
1123
1124            for (i = 0; i < l; i++)
1125                alpha[i] *= y[i] / r;
1126
1127            si.rho /= r;
1128            si.obj /= (r * r);
1129            si.upper_bound_p = 1 / r;
1130            si.upper_bound_n = 1 / r;
1131        }
1132
1133        private static void solve_one_class(Problem prob, Parameter param,
1134                        double[] alpha, Solver.SolutionInfo si)
1135        {
1136            int l = prob.Count;
1137            double[] zeros = new double[l];
1138            sbyte[] ones = new sbyte[l];
1139            int i;
1140
1141            int n = (int)(param.Nu * prob.Count); // # of alpha's at upper bound
1142
1143            for (i = 0; i < n; i++)
1144                alpha[i] = 1;
1145            if (n < prob.Count)
1146                alpha[n] = param.Nu * prob.Count - n;
1147            for (i = n + 1; i < l; i++)
1148                alpha[i] = 0;
1149
1150            for (i = 0; i < l; i++)
1151            {
1152                zeros[i] = 0;
1153                ones[i] = 1;
1154            }
1155
1156            Solver s = new Solver();
1157            s.Solve(l, new ONE_CLASS_Q(prob, param), zeros, ones, alpha, 1.0, 1.0, param.EPS, si, param.Shrinking);
1158        }
1159
1160        private static void solve_epsilon_svr(Problem prob, Parameter param, double[] alpha, Solver.SolutionInfo si)
1161        {
1162            int l = prob.Count;
1163            double[] alpha2 = new double[2 * l];
1164            double[] linear_term = new double[2 * l];
1165            sbyte[] y = new sbyte[2 * l];
1166            int i;
1167
1168            for (i = 0; i < l; i++)
1169            {
1170                alpha2[i] = 0;
1171                linear_term[i] = param.P - prob.Y[i];
1172                y[i] = 1;
1173
1174                alpha2[i + l] = 0;
1175                linear_term[i + l] = param.P + prob.Y[i];
1176                y[i + l] = -1;
1177            }
1178
1179            Solver s = new Solver();
1180            s.Solve(2 * l, new SVR_Q(prob, param), linear_term, y, alpha2, param.C, param.C, param.EPS, si, param.Shrinking);
1181
1182            double sum_alpha = 0;
1183            for (i = 0; i < l; i++)
1184            {
1185                alpha[i] = alpha2[i] - alpha2[i + l];
1186                sum_alpha += Math.Abs(alpha[i]);
1187            }
1188            Procedures.info("nu = " + sum_alpha / (param.C * l) + "\n");
1189        }
1190
1191        private static void solve_nu_svr(Problem prob, Parameter param,
1192                        double[] alpha, Solver.SolutionInfo si)
1193        {
1194            int l = prob.Count;
1195            double C = param.C;
1196            double[] alpha2 = new double[2 * l];
1197            double[] linear_term = new double[2 * l];
1198            sbyte[] y = new sbyte[2 * l];
1199            int i;
1200
1201            double sum = C * param.Nu * l / 2;
1202            for (i = 0; i < l; i++)
1203            {
1204                alpha2[i] = alpha2[i + l] = Math.Min(sum, C);
1205                sum -= alpha2[i];
1206
1207                linear_term[i] = -prob.Y[i];
1208                y[i] = 1;
1209
1210                linear_term[i + l] = prob.Y[i];
1211                y[i + l] = -1;
1212            }
1213
1214            Solver_NU s = new Solver_NU();
1215            s.Solve(2 * l, new SVR_Q(prob, param), linear_term, y, alpha2, C, C, param.EPS, si, param.Shrinking);
1216
1217            Procedures.info("epsilon = " + (-si.r) + "\n");
1218
1219            for (i = 0; i < l; i++)
1220                alpha[i] = alpha2[i] - alpha2[i + l];
1221        }
1222
1223        //
1224        // decision_function
1225        //
1226        internal class decision_function
1227        {
1228            public double[] alpha;
1229            public double rho;
1230        };
1231
1232        static decision_function svm_train_one(Problem prob, Parameter param, double Cp, double Cn)
1233        {
1234            double[] alpha = new double[prob.Count];
1235            Solver.SolutionInfo si = new Solver.SolutionInfo();
1236            switch (param.SvmType)
1237            {
1238                case SvmType.C_SVC:
1239                    solve_c_svc(prob, param, alpha, si, Cp, Cn);
1240                    break;
1241                case SvmType.NU_SVC:
1242                    solve_nu_svc(prob, param, alpha, si);
1243                    break;
1244                case SvmType.ONE_CLASS:
1245                    solve_one_class(prob, param, alpha, si);
1246                    break;
1247                case SvmType.EPSILON_SVR:
1248                    solve_epsilon_svr(prob, param, alpha, si);
1249                    break;
1250                case SvmType.NU_SVR:
1251                    solve_nu_svr(prob, param, alpha, si);
1252                    break;
1253            }
1254
1255            Procedures.info("obj = " + si.obj + ", rho = " + si.rho + "\n");
1256
1257            // output SVs
1258
1259            int nSV = 0;
1260            int nBSV = 0;
1261            for (int i = 0; i < prob.Count; i++)
1262            {
1263                if (Math.Abs(alpha[i]) > 0)
1264                {
1265                    ++nSV;
1266                    if (prob.Y[i] > 0)
1267                    {
1268                        if (Math.Abs(alpha[i]) >= si.upper_bound_p)
1269                            ++nBSV;
1270                    }
1271                    else
1272                    {
1273                        if (Math.Abs(alpha[i]) >= si.upper_bound_n)
1274                            ++nBSV;
1275                    }
1276                }
1277            }
1278
1279            Procedures.info("nSV = " + nSV + ", nBSV = " + nBSV + "\n");
1280
1281            decision_function f = new decision_function();
1282            f.alpha = alpha;
1283            f.rho = si.rho;
1284            return f;
1285        }
1286
1287        // Platt's binary SVM Probablistic Output: an improvement from Lin et al.
1288        private static void sigmoid_train(int l, double[] dec_values, double[] labels,
1289                      double[] probAB)
1290        {
1291            double A, B;
1292            double prior1 = 0, prior0 = 0;
1293            int i;
1294
1295            for (i = 0; i < l; i++)
1296                if (labels[i] > 0) prior1 += 1;
1297                else prior0 += 1;
1298
1299            int Max_iter = 100; // Maximal number of iterations
1300            double Min_step = 1e-10;  // Minimal step taken in line search
1301            double sigma = 1e-12; // For numerically strict PD of Hessian
1302            double eps = 1e-5;
1303            double hiTarget = (prior1 + 1.0) / (prior1 + 2.0);
1304            double loTarget = 1 / (prior0 + 2.0);
1305            double[] t = new double[l];
1306            double fApB, p, q, h11, h22, h21, g1, g2, det, dA, dB, gd, stepsize;
1307            double newA, newB, newf, d1, d2;
1308            int iter;
1309
1310            // Initial Point and Initial Fun Value
1311            A = 0.0; B = Math.Log((prior0 + 1.0) / (prior1 + 1.0));
1312            double fval = 0.0;
1313
1314            for (i = 0; i < l; i++)
1315            {
1316                if (labels[i] > 0) t[i] = hiTarget;
1317                else t[i] = loTarget;
1318                fApB = dec_values[i] * A + B;
1319                if (fApB >= 0)
1320                    fval += t[i] * fApB + Math.Log(1 + Math.Exp(-fApB));
1321                else
1322                    fval += (t[i] - 1) * fApB + Math.Log(1 + Math.Exp(fApB));
1323            }
1324            for (iter = 0; iter < Max_iter; iter++)
1325            {
1326                // Update Gradient and Hessian (use H' = H + sigma I)
1327                h11 = sigma; // numerically ensures strict PD
1328                h22 = sigma;
1329                h21 = 0.0; g1 = 0.0; g2 = 0.0;
1330                for (i = 0; i < l; i++)
1331                {
1332                    fApB = dec_values[i] * A + B;
1333                    if (fApB >= 0)
1334                    {
1335                        p = Math.Exp(-fApB) / (1.0 + Math.Exp(-fApB));
1336                        q = 1.0 / (1.0 + Math.Exp(-fApB));
1337                    }
1338                    else
1339                    {
1340                        p = 1.0 / (1.0 + Math.Exp(fApB));
1341                        q = Math.Exp(fApB) / (1.0 + Math.Exp(fApB));
1342                    }
1343                    d2 = p * q;
1344                    h11 += dec_values[i] * dec_values[i] * d2;
1345                    h22 += d2;
1346                    h21 += dec_values[i] * d2;
1347                    d1 = t[i] - p;
1348                    g1 += dec_values[i] * d1;
1349                    g2 += d1;
1350                }
1351
1352                // Stopping Criteria
1353                if (Math.Abs(g1) < eps && Math.Abs(g2) < eps)
1354                    break;
1355
1356                // Finding Newton direction: -inv(H') * g
1357                det = h11 * h22 - h21 * h21;
1358                dA = -(h22 * g1 - h21 * g2) / det;
1359                dB = -(-h21 * g1 + h11 * g2) / det;
1360                gd = g1 * dA + g2 * dB;
1361
1362
1363                stepsize = 1;   // Line Search
1364                while (stepsize >= Min_step)
1365                {
1366                    newA = A + stepsize * dA;
1367                    newB = B + stepsize * dB;
1368
1369                    // New function value
1370                    newf = 0.0;
1371                    for (i = 0; i < l; i++)
1372                    {
1373                        fApB = dec_values[i] * newA + newB;
1374                        if (fApB >= 0)
1375                            newf += t[i] * fApB + Math.Log(1 + Math.Exp(-fApB));
1376                        else
1377                            newf += (t[i] - 1) * fApB + Math.Log(1 + Math.Exp(fApB));
1378                    }
1379                    // Check sufficient decrease
1380                    if (newf < fval + 0.0001 * stepsize * gd)
1381                    {
1382                        A = newA; B = newB; fval = newf;
1383                        break;
1384                    }
1385                    else
1386                        stepsize = stepsize / 2.0;
1387                }
1388
1389                if (stepsize < Min_step)
1390                {
1391                    Procedures.info("Line search fails in two-class probability estimates\n");
1392                    break;
1393                }
1394            }
1395
1396            if (iter >= Max_iter)
1397                Procedures.info("Reaching Maximal iterations in two-class probability estimates\n");
1398            probAB[0] = A; probAB[1] = B;
1399        }
1400
1401        private static double sigmoid_predict(double decision_value, double A, double B)
1402        {
1403            double fApB = decision_value * A + B;
1404            if (fApB >= 0)
1405                return Math.Exp(-fApB) / (1.0 + Math.Exp(-fApB));
1406            else
1407                return 1.0 / (1 + Math.Exp(fApB));
1408        }
1409
1410        // Method 2 from the multiclass_prob paper by Wu, Lin, and Weng
1411        private static void multiclass_probability(int k, double[,] r, double[] p)
1412        {
1413            int t, j;
1414            int iter = 0, Max_iter = Math.Max(100, k);
1415            double[,] Q = new double[k,k];
1416            double[] Qp = new double[k];
1417            double pQp, eps = 0.005 / k;
1418
1419            for (t = 0; t < k; t++)
1420            {
1421                p[t] = 1.0 / k;  // Valid if k = 1
1422                Q[t,t] = 0;
1423                for (j = 0; j < t; j++)
1424                {
1425                    Q[t,t] += r[j,t] * r[j,t];
1426                    Q[t,j] = Q[j,t];
1427                }
1428                for (j = t + 1; j < k; j++)
1429                {
1430                    Q[t,t] += r[j,t] * r[j,t];
1431                    Q[t,j] = -r[j,t] * r[t,j];
1432                }
1433            }
1434            for (iter = 0; iter < Max_iter; iter++)
1435            {
1436                // stopping condition, recalculate QP,pQP for numerical accuracy
1437                pQp = 0;
1438                for (t = 0; t < k; t++)
1439                {
1440                    Qp[t] = 0;
1441                    for (j = 0; j < k; j++)
1442                        Qp[t] += Q[t,j] * p[j];
1443                    pQp += p[t] * Qp[t];
1444                }
1445                double Max_error = 0;
1446                for (t = 0; t < k; t++)
1447                {
1448                    double error = Math.Abs(Qp[t] - pQp);
1449                    if (error > Max_error)
1450                        Max_error = error;
1451                }
1452                if (Max_error < eps) break;
1453
1454                for (t = 0; t < k; t++)
1455                {
1456                    double diff = (-Qp[t] + pQp) / Q[t,t];
1457                    p[t] += diff;
1458                    pQp = (pQp + diff * (diff * Q[t,t] + 2 * Qp[t])) / (1 + diff) / (1 + diff);
1459                    for (j = 0; j < k; j++)
1460                    {
1461                        Qp[j] = (Qp[j] + diff * Q[t,j]) / (1 + diff);
1462                        p[j] /= (1 + diff);
1463                    }
1464                }
1465            }
1466            if (iter >= Max_iter)
1467                Procedures.info("Exceeds Max_iter in multiclass_prob\n");
1468        }
1469
1470        // Cross-validation decision values for probability estimates
1471        private static void svm_binary_svc_probability(Problem prob, Parameter param, double Cp, double Cn, double[] probAB)
1472        {
1473            int i;
1474            int nr_fold = 5;
1475            int[] perm = new int[prob.Count];
1476            double[] dec_values = new double[prob.Count];
1477
1478            // random shuffle
1479            Random rand = new Random();
1480            for (i = 0; i < prob.Count; i++) perm[i] = i;
1481            for (i = 0; i < prob.Count; i++)
1482            {
1483                int j = i + (int)(rand.NextDouble() * (prob.Count - i));
1484                do { int _ = perm[i]; perm[i] = perm[j]; perm[j] = _; } while (false);
1485            }
1486            for (i = 0; i < nr_fold; i++)
1487            {
1488                int begin = i * prob.Count / nr_fold;
1489                int end = (i + 1) * prob.Count / nr_fold;
1490                int j, k;
1491                Problem subprob = new Problem();
1492
1493                subprob.Count = prob.Count - (end - begin);
1494                subprob.X = new Node[subprob.Count][];
1495                subprob.Y = new double[subprob.Count];
1496
1497                k = 0;
1498                for (j = 0; j < begin; j++)
1499                {
1500                    subprob.X[k] = prob.X[perm[j]];
1501                    subprob.Y[k] = prob.Y[perm[j]];
1502                    ++k;
1503                }
1504                for (j = end; j < prob.Count; j++)
1505                {
1506                    subprob.X[k] = prob.X[perm[j]];
1507                    subprob.Y[k] = prob.Y[perm[j]];
1508                    ++k;
1509                }
1510                int p_count = 0, n_count = 0;
1511                for (j = 0; j < k; j++)
1512                    if (subprob.Y[j] > 0)
1513                        p_count++;
1514                    else
1515                        n_count++;
1516
1517                if (p_count == 0 && n_count == 0)
1518                    for (j = begin; j < end; j++)
1519                        dec_values[perm[j]] = 0;
1520                else if (p_count > 0 && n_count == 0)
1521                    for (j = begin; j < end; j++)
1522                        dec_values[perm[j]] = 1;
1523                else if (p_count == 0 && n_count > 0)
1524                    for (j = begin; j < end; j++)
1525                        dec_values[perm[j]] = -1;
1526                else
1527                {
1528                    Parameter subparam = (Parameter)param.Clone();
1529                    subparam.Probability = false;
1530                    subparam.C = 1.0;
1531                    subparam.Weights[1] = Cp;
1532                    subparam.Weights[-1] = Cn;
1533                    Model submodel = svm_train(subprob, subparam);
1534                    for (j = begin; j < end; j++)
1535                    {
1536                        double[] dec_value = new double[1];
1537                        svm_predict_values(submodel, prob.X[perm[j]], dec_value);
1538                        dec_values[perm[j]] = dec_value[0];
1539                        // ensure +1 -1 order; reason not using CV subroutine
1540                        dec_values[perm[j]] *= submodel.ClassLabels[0];
1541                    }
1542                }
1543            }
1544            sigmoid_train(prob.Count, dec_values, prob.Y, probAB);
1545        }
1546
1547        // Return parameter of a Laplace distribution
1548        private static double svm_svr_probability(Problem prob, Parameter param)
1549        {
1550            int i;
1551            int nr_fold = 5;
1552            double[] ymv = new double[prob.Count];
1553            double mae = 0;
1554
1555            Parameter newparam = (Parameter)param.Clone();
1556            newparam.Probability = false;
1557            svm_cross_validation(prob, newparam, nr_fold, ymv, true);
1558            for (i = 0; i < prob.Count; i++)
1559            {
1560                ymv[i] = prob.Y[i] - ymv[i];
1561                mae += Math.Abs(ymv[i]);
1562            }
1563            mae /= prob.Count;
1564            double std = Math.Sqrt(2 * mae * mae);
1565            int count = 0;
1566            mae = 0;
1567            for (i = 0; i < prob.Count; i++)
1568                if (Math.Abs(ymv[i]) > 5 * std)
1569                    count = count + 1;
1570                else
1571                    mae += Math.Abs(ymv[i]);
1572            mae /= (prob.Count - count);
1573            Procedures.info("Prob. model for test data: target value = predicted value + z,\nz: Laplace distribution e^(-|z|/sigma)/(2sigma),sigma=" + mae + "\n");
1574            return mae;
1575        }
1576
1577        // label: label name, start: begin of each class, count: #data of classes, perm: indices to the original data
1578        // perm, length l, must be allocated before calling this subroutine
1579        private static void svm_group_classes(Problem prob, int[] nr_class_ret, int[][] label_ret, int[][] start_ret, int[][] count_ret, int[] perm)
1580        {
1581            int l = prob.Count;
1582            int Max_nr_class = 16;
1583            int nr_class = 0;
1584            int[] label = new int[Max_nr_class];
1585            int[] count = new int[Max_nr_class];
1586            int[] data_label = new int[l];
1587            int i;
1588
1589            for (i = 0; i < l; i++)
1590            {
1591                int this_label = (int)(prob.Y[i]);
1592                int j;
1593                for (j = 0; j < nr_class; j++)
1594                {
1595                    if (this_label == label[j])
1596                    {
1597                        ++count[j];
1598                        break;
1599                    }
1600                }
1601                data_label[i] = j;
1602                if (j == nr_class)
1603                {
1604                    if (nr_class == Max_nr_class)
1605                    {
1606                        Max_nr_class *= 2;
1607                        int[] new_data = new int[Max_nr_class];
1608                        Array.Copy(label, 0, new_data, 0, label.Length);
1609                        label = new_data;
1610                        new_data = new int[Max_nr_class];
1611                        Array.Copy(count, 0, new_data, 0, count.Length);
1612                        count = new_data;
1613                    }
1614                    label[nr_class] = this_label;
1615                    count[nr_class] = 1;
1616                    ++nr_class;
1617                }
1618            }
1619
1620            int[] start = new int[nr_class];
1621            start[0] = 0;
1622            for (i = 1; i < nr_class; i++)
1623                start[i] = start[i - 1] + count[i - 1];
1624            for (i = 0; i < l; i++)
1625            {
1626                perm[start[data_label[i]]] = i;
1627                ++start[data_label[i]];
1628            }
1629            start[0] = 0;
1630            for (i = 1; i < nr_class; i++)
1631                start[i] = start[i - 1] + count[i - 1];
1632
1633            nr_class_ret[0] = nr_class;
1634            label_ret[0] = label;
1635            start_ret[0] = start;
1636            count_ret[0] = count;
1637        }
1638
1639        //
1640        // Interface functions
1641        //
1642        public static Model svm_train(Problem prob, Parameter param)
1643        {
1644            Model model = new Model();
1645            model.Parameter = param;
1646
1647            if (param.SvmType == SvmType.ONE_CLASS ||
1648               param.SvmType == SvmType.EPSILON_SVR ||
1649               param.SvmType == SvmType.NU_SVR)
1650            {
1651                // regression or one-class-svm
1652                model.NumberOfClasses = 2;               
1653                model.ClassLabels = null;
1654                model.NumberOfSVPerClass = null;
1655                model.PairwiseProbabilityA = null; model.PairwiseProbabilityB = null;
1656                model.SupportVectorCoefficients = new double[1][];
1657               
1658                if (param.Probability &&
1659                   (param.SvmType == SvmType.EPSILON_SVR ||
1660                    param.SvmType == SvmType.NU_SVR))
1661                {
1662                    model.PairwiseProbabilityA = new double[1];
1663                    model.PairwiseProbabilityA[0] = svm_svr_probability(prob, param);
1664                }
1665
1666                decision_function f = svm_train_one(prob, param, 0, 0);
1667                model.Rho = new double[1];
1668                model.Rho[0] = f.rho;
1669
1670                int nSV = 0;
1671                int i;
1672                for (i = 0; i < prob.Count; i++)
1673                    if (Math.Abs(f.alpha[i]) > 0) ++nSV;               
1674                model.SupportVectorCount = nSV;
1675                model.SupportVectorIndizes = new int[nSV];
1676                model.SupportVectors = new Node[nSV][];
1677                model.SupportVectorCoefficients[0] = new double[nSV];
1678               
1679                int j = 0;
1680                for (i = 0; i < prob.Count; i++)
1681                    if (Math.Abs(f.alpha[i]) > 0)
1682                    {
1683                        model.SupportVectors[j] = prob.X[i];
1684                        model.SupportVectorIndizes[j] = i;
1685                        model.SupportVectorCoefficients[0][j] = f.alpha[i];
1686                       
1687                        ++j;
1688                    }
1689            }
1690            else
1691            {
1692                // classification
1693                int l = prob.Count;
1694                int[] tmp_nr_class = new int[1];
1695                int[][] tmp_label = new int[1][];
1696                int[][] tmp_start = new int[1][];
1697                int[][] tmp_count = new int[1][];
1698                int[] perm = new int[l];
1699
1700                // group training data of the same class
1701                svm_group_classes(prob, tmp_nr_class, tmp_label, tmp_start, tmp_count, perm);
1702                int nr_class = tmp_nr_class[0];
1703                int[] label = tmp_label[0];
1704                int[] start = tmp_start[0];
1705                int[] count = tmp_count[0];
1706                Node[][] x = new Node[l][];
1707                int i;
1708                for (i = 0; i < l; i++)
1709                    x[i] = prob.X[perm[i]];
1710
1711                // calculate weighted C
1712
1713                double[] weighted_C = new double[nr_class];
1714                for (i = 0; i < nr_class; i++)
1715                    weighted_C[i] = param.C;
1716                foreach (int weightedLabel in param.Weights.Keys)
1717                {
1718                    int index = Array.IndexOf<int>(label, weightedLabel);
1719                    if (index < 0)
1720                        Console.Error.WriteLine("warning: class label " + weightedLabel + " specified in weight is not found");
1721                    else weighted_C[index] *= param.Weights[weightedLabel];
1722                }
1723
1724                // train k*(k-1)/2 models
1725
1726                bool[] nonzero = new bool[l];
1727                for (i = 0; i < l; i++)
1728                    nonzero[i] = false;
1729                decision_function[] f = new decision_function[nr_class * (nr_class - 1) / 2];
1730
1731                double[] probA = null, probB = null;
1732                if (param.Probability)
1733                {
1734                    probA = new double[nr_class * (nr_class - 1) / 2];
1735                    probB = new double[nr_class * (nr_class - 1) / 2];
1736                }
1737
1738                int p = 0;
1739                for (i = 0; i < nr_class; i++)
1740                    for (int j = i + 1; j < nr_class; j++)
1741                    {
1742                        Problem sub_prob = new Problem();
1743                        int si = start[i], sj = start[j];
1744                        int ci = count[i], cj = count[j];
1745                        sub_prob.Count = ci + cj;
1746                        sub_prob.X = new Node[sub_prob.Count][];
1747                        sub_prob.Y = new double[sub_prob.Count];
1748                        int k;
1749                        for (k = 0; k < ci; k++)
1750                        {
1751                            sub_prob.X[k] = x[si + k];
1752                            sub_prob.Y[k] = +1;
1753                        }
1754                        for (k = 0; k < cj; k++)
1755                        {
1756                            sub_prob.X[ci + k] = x[sj + k];
1757                            sub_prob.Y[ci + k] = -1;
1758                        }
1759
1760                        if (param.Probability)
1761                        {
1762                            double[] probAB = new double[2];
1763                            svm_binary_svc_probability(sub_prob, param, weighted_C[i], weighted_C[j], probAB);
1764                            probA[p] = probAB[0];
1765                            probB[p] = probAB[1];
1766                        }
1767
1768                        f[p] = svm_train_one(sub_prob, param, weighted_C[i], weighted_C[j]);
1769                        for (k = 0; k < ci; k++)
1770                            if (!nonzero[si + k] && Math.Abs(f[p].alpha[k]) > 0)
1771                                nonzero[si + k] = true;
1772                        for (k = 0; k < cj; k++)
1773                            if (!nonzero[sj + k] && Math.Abs(f[p].alpha[ci + k]) > 0)
1774                                nonzero[sj + k] = true;
1775                        ++p;
1776                    }
1777
1778                // build output
1779
1780                model.NumberOfClasses = nr_class;
1781
1782                model.ClassLabels = new int[nr_class];
1783                for (i = 0; i < nr_class; i++)
1784                    model.ClassLabels[i] = label[i];
1785
1786                model.Rho = new double[nr_class * (nr_class - 1) / 2];
1787                for (i = 0; i < nr_class * (nr_class - 1) / 2; i++)
1788                    model.Rho[i] = f[i].rho;
1789
1790                if (param.Probability)
1791                {
1792                    model.PairwiseProbabilityA = new double[nr_class * (nr_class - 1) / 2];
1793                    model.PairwiseProbabilityB = new double[nr_class * (nr_class - 1) / 2];
1794                    for (i = 0; i < nr_class * (nr_class - 1) / 2; i++)
1795                    {
1796                        model.PairwiseProbabilityA[i] = probA[i];
1797                        model.PairwiseProbabilityB[i] = probB[i];
1798                    }
1799                }
1800                else
1801                {
1802                    model.PairwiseProbabilityA = null;
1803                    model.PairwiseProbabilityB = null;
1804                }
1805
1806                int nnz = 0;
1807                int[] nz_count = new int[nr_class];
1808                model.NumberOfSVPerClass = new int[nr_class];
1809                for (i = 0; i < nr_class; i++)
1810                {
1811                    int nSV = 0;
1812                    for (int j = 0; j < count[i]; j++)
1813                        if (nonzero[start[i] + j])
1814                        {
1815                            ++nSV;
1816                            ++nnz;
1817                        }
1818                    model.NumberOfSVPerClass[i] = nSV;
1819                    nz_count[i] = nSV;
1820                }
1821
1822                Procedures.info("Total nSV = " + nnz + "\n");
1823
1824                model.SupportVectorCount = nnz;
1825                model.SupportVectors = new Node[nnz][];
1826                model.SupportVectorIndizes = new int[nnz];
1827                p = 0;
1828                for (i = 0; i < l; i++) {
1829                  if (nonzero[i]) {
1830                    model.SupportVectors[p] = x[i];
1831                    model.SupportVectorIndizes[p] = i;
1832                    p++;
1833                  }
1834                }
1835
1836                int[] nz_start = new int[nr_class];
1837                nz_start[0] = 0;
1838                for (i = 1; i < nr_class; i++)
1839                    nz_start[i] = nz_start[i - 1] + nz_count[i - 1];
1840
1841                model.SupportVectorCoefficients = new double[nr_class - 1][];
1842                for (i = 0; i < nr_class - 1; i++)
1843                    model.SupportVectorCoefficients[i] = new double[nnz];
1844
1845                p = 0;
1846                for (i = 0; i < nr_class; i++)
1847                    for (int j = i + 1; j < nr_class; j++)
1848                    {
1849                        // classifier (i,j): coefficients with
1850                        // i are in sv_coef[j-1][nz_start[i]...],
1851                        // j are in sv_coef[i][nz_start[j]...]
1852
1853                        int si = start[i];
1854                        int sj = start[j];
1855                        int ci = count[i];
1856                        int cj = count[j];
1857
1858                        int q = nz_start[i];
1859                        int k;
1860                        for (k = 0; k < ci; k++)
1861                            if (nonzero[si + k])
1862                                model.SupportVectorCoefficients[j - 1][q++] = f[p].alpha[k];
1863                        q = nz_start[j];
1864                        for (k = 0; k < cj; k++)
1865                            if (nonzero[sj + k])
1866                                model.SupportVectorCoefficients[i][q++] = f[p].alpha[ci + k];
1867                        ++p;
1868                    }
1869            }
1870            return model;
1871        }
1872
1873        // Stratified cross validation
1874        public static void svm_cross_validation(Problem prob, Parameter param, int nr_fold, double[] target, bool shuffleTraining)
1875        {
1876            Random rand = new Random();
1877            int i;
1878            int[] fold_start = new int[nr_fold + 1];
1879            int l = prob.Count;
1880            int[] perm = new int[l];
1881
1882            // stratified cv may not give leave-one-out rate
1883            // Each class to l folds -> some folds may have zero elements
1884            if ((param.SvmType == SvmType.C_SVC ||
1885                param.SvmType == SvmType.NU_SVC) && nr_fold < l)
1886            {
1887                int[] tmp_nr_class = new int[1];
1888                int[][] tmp_label = new int[1][];
1889                int[][] tmp_start = new int[1][];
1890                int[][] tmp_count = new int[1][];
1891
1892                svm_group_classes(prob, tmp_nr_class, tmp_label, tmp_start, tmp_count, perm);
1893
1894                int nr_class = tmp_nr_class[0];
1895                int[] label = tmp_label[0];
1896                int[] start = tmp_start[0];
1897                int[] count = tmp_count[0];
1898
1899                // random shuffle and then data grouped by fold using the array perm
1900                int[] fold_count = new int[nr_fold];
1901                int c;
1902                int[] index = new int[l];
1903                for (i = 0; i < l; i++)
1904                    index[i] = perm[i];
1905                for (c = 0; c < nr_class; c++)
1906                    for (i = 0; i < count[c]; i++)
1907                    {
1908                        int j = i + (int)(rand.NextDouble() * (count[c] - i));
1909                        do { int _ = index[start[c] + j]; index[start[c] + j] = index[start[c] + i]; index[start[c] + i] = _; } while (false);
1910                    }
1911                for (i = 0; i < nr_fold; i++)
1912                {
1913                    fold_count[i] = 0;
1914                    for (c = 0; c < nr_class; c++)
1915                        fold_count[i] += (i + 1) * count[c] / nr_fold - i * count[c] / nr_fold;
1916                }
1917                fold_start[0] = 0;
1918                for (i = 1; i <= nr_fold; i++)
1919                    fold_start[i] = fold_start[i - 1] + fold_count[i - 1];
1920                for (c = 0; c < nr_class; c++)
1921                    for (i = 0; i < nr_fold; i++)
1922                    {
1923                        int begin = start[c] + i * count[c] / nr_fold;
1924                        int end = start[c] + (i + 1) * count[c] / nr_fold;
1925                        for (int j = begin; j < end; j++)
1926                        {
1927                            perm[fold_start[i]] = index[j];
1928                            fold_start[i]++;
1929                        }
1930                    }
1931                fold_start[0] = 0;
1932                for (i = 1; i <= nr_fold; i++)
1933                    fold_start[i] = fold_start[i - 1] + fold_count[i - 1];
1934            }
1935            else
1936            {
1937                for (i = 0; i < l; i++) perm[i] = i;
1938                if (shuffleTraining) {
1939                  for (i = 0; i < l; i++) {
1940                    int j = i + (int)(rand.NextDouble() * (l - i));
1941                    do { int _ = perm[i]; perm[i] = perm[j]; perm[j] = _; } while (false);
1942                  }
1943                }
1944                for (i = 0; i <= nr_fold; i++)
1945                    fold_start[i] = i * l / nr_fold;
1946            }
1947
1948            for (i = 0; i < nr_fold; i++)
1949            {
1950                int begin = fold_start[i];
1951                int end = fold_start[i + 1];
1952                int j, k;
1953                Problem subprob = new Problem();
1954
1955                subprob.Count = l - (end - begin);
1956                subprob.X = new Node[subprob.Count][];
1957                subprob.Y = new double[subprob.Count];
1958
1959                k = 0;
1960                for (j = 0; j < begin; j++)
1961                {
1962                    subprob.X[k] = prob.X[perm[j]];
1963                    subprob.Y[k] = prob.Y[perm[j]];
1964                    ++k;
1965                }
1966                for (j = end; j < l; j++)
1967                {
1968                    subprob.X[k] = prob.X[perm[j]];
1969                    subprob.Y[k] = prob.Y[perm[j]];
1970                    ++k;
1971                }
1972                Model submodel = svm_train(subprob, param);
1973                if (param.Probability &&
1974                   (param.SvmType == SvmType.C_SVC ||
1975                    param.SvmType == SvmType.NU_SVC))
1976                {
1977                    double[] prob_estimates = new double[svm_get_nr_class(submodel)];
1978                    for (j = begin; j < end; j++)
1979                        target[perm[j]] = svm_predict_probability(submodel, prob.X[perm[j]], prob_estimates);
1980                }
1981                else
1982                    for (j = begin; j < end; j++)
1983                        target[perm[j]] = svm_predict(submodel, prob.X[perm[j]]);
1984            }
1985        }
1986
1987        public static SvmType svm_get_svm_type(Model model)
1988        {
1989            return model.Parameter.SvmType;
1990        }
1991
1992        public static int svm_get_nr_class(Model model)
1993        {
1994            return model.NumberOfClasses;
1995        }
1996
1997        public static void svm_get_labels(Model model, int[] label)
1998        {
1999            if (model.ClassLabels != null)
2000                for (int i = 0; i < model.NumberOfClasses; i++)
2001                    label[i] = model.ClassLabels[i];
2002        }
2003
2004        public static double svm_get_svr_probability(Model model)
2005        {
2006            if ((model.Parameter.SvmType == SvmType.EPSILON_SVR || model.Parameter.SvmType == SvmType.NU_SVR) &&
2007                model.PairwiseProbabilityA != null)
2008                return model.PairwiseProbabilityA[0];
2009            else
2010            {
2011                Console.Error.WriteLine("Model doesn't contain information for SVR probability inference");
2012                return 0;
2013            }
2014        }
2015
2016        public static void svm_predict_values(Model model, Node[] x, double[] dec_values)
2017        {
2018            if (model.Parameter.SvmType == SvmType.ONE_CLASS ||
2019               model.Parameter.SvmType == SvmType.EPSILON_SVR ||
2020               model.Parameter.SvmType == SvmType.NU_SVR)
2021            {
2022                double[] sv_coef = model.SupportVectorCoefficients[0];
2023                double sum = 0;
2024                for (int i = 0; i < model.SupportVectorCount; i++)
2025                    sum += sv_coef[i] * Kernel.KernelFunction(x, model.SupportVectors[i], model.Parameter);
2026                sum -= model.Rho[0];
2027                dec_values[0] = sum;
2028            }
2029            else
2030            {
2031                int i;
2032                int nr_class = model.NumberOfClasses;
2033                int l = model.SupportVectorCount;
2034
2035                double[] kvalue = new double[l];
2036                for (i = 0; i < l; i++)
2037                    kvalue[i] = Kernel.KernelFunction(x, model.SupportVectors[i], model.Parameter);
2038
2039                int[] start = new int[nr_class];
2040                start[0] = 0;
2041                for (i = 1; i < nr_class; i++)
2042                    start[i] = start[i - 1] + model.NumberOfSVPerClass[i - 1];
2043
2044                int p = 0;
2045                for (i = 0; i < nr_class; i++)
2046                    for (int j = i + 1; j < nr_class; j++)
2047                    {
2048                        double sum = 0;
2049                        int si = start[i];
2050                        int sj = start[j];
2051                        int ci = model.NumberOfSVPerClass[i];
2052                        int cj = model.NumberOfSVPerClass[j];
2053
2054                        int k;
2055                        double[] coef1 = model.SupportVectorCoefficients[j - 1];
2056                        double[] coef2 = model.SupportVectorCoefficients[i];
2057                        for (k = 0; k < ci; k++)
2058                            sum += coef1[si + k] * kvalue[si + k];
2059                        for (k = 0; k < cj; k++)
2060                            sum += coef2[sj + k] * kvalue[sj + k];
2061                        sum -= model.Rho[p];
2062                        dec_values[p] = sum;
2063                        p++;
2064                    }
2065            }
2066        }
2067
2068        public static double svm_predict(Model model, Node[] x)
2069        {
2070            if (model.Parameter.SvmType == SvmType.ONE_CLASS ||
2071               model.Parameter.SvmType == SvmType.EPSILON_SVR ||
2072               model.Parameter.SvmType == SvmType.NU_SVR)
2073            {
2074                double[] res = new double[1];
2075                svm_predict_values(model, x, res);
2076
2077                if (model.Parameter.SvmType == SvmType.ONE_CLASS)
2078                    return (res[0] > 0) ? 1 : -1;
2079                else
2080                    return res[0];
2081            }
2082            else
2083            {
2084                int i;
2085                int nr_class = model.NumberOfClasses;
2086                double[] dec_values = new double[nr_class * (nr_class - 1) / 2];
2087                svm_predict_values(model, x, dec_values);
2088
2089                int[] vote = new int[nr_class];
2090                for (i = 0; i < nr_class; i++)
2091                    vote[i] = 0;
2092                int pos = 0;
2093                for (i = 0; i < nr_class; i++)
2094                    for (int j = i + 1; j < nr_class; j++)
2095                    {
2096                        if (dec_values[pos++] > 0)
2097                            ++vote[i];
2098                        else
2099                            ++vote[j];
2100                    }
2101
2102                int vote_Max_idx = 0;
2103                for (i = 1; i < nr_class; i++)
2104                    if (vote[i] > vote[vote_Max_idx])
2105                        vote_Max_idx = i;
2106                return model.ClassLabels[vote_Max_idx];
2107            }
2108        }
2109
2110        public static double svm_predict_probability(Model model, Node[] x, double[] prob_estimates)
2111        {
2112            if ((model.Parameter.SvmType == SvmType.C_SVC || model.Parameter.SvmType == SvmType.NU_SVC) &&
2113                model.PairwiseProbabilityA != null && model.PairwiseProbabilityB != null)
2114            {
2115                int i;
2116                int nr_class = model.NumberOfClasses;
2117                double[] dec_values = new double[nr_class * (nr_class - 1) / 2];
2118                svm_predict_values(model, x, dec_values);
2119
2120                double Min_prob = 1e-7;
2121                double[,] pairwise_prob = new double[nr_class, nr_class];
2122
2123                int k = 0;
2124                for (i = 0; i < nr_class; i++)
2125                {
2126                    for (int j = i + 1; j < nr_class; j++)
2127                    {
2128                        pairwise_prob[i, j] = Math.Min(Math.Max(sigmoid_predict(dec_values[k], model.PairwiseProbabilityA[k], model.PairwiseProbabilityB[k]), Min_prob), 1 - Min_prob);
2129                        pairwise_prob[j, i] = 1 - pairwise_prob[i, j];
2130                        k++;
2131                    }
2132                }
2133                multiclass_probability(nr_class, pairwise_prob, prob_estimates);
2134
2135                int prob_Max_idx = 0;
2136                for (i = 1; i < nr_class; i++)
2137                    if (prob_estimates[i] > prob_estimates[prob_Max_idx])
2138                        prob_Max_idx = i;
2139                return model.ClassLabels[prob_Max_idx];
2140            }
2141            else
2142                return svm_predict(model, x);
2143        }
2144
2145        public static string svm_check_parameter(Problem prob, Parameter param)
2146        {
2147            // svm_type
2148
2149            SvmType svm_type = param.SvmType;
2150
2151            // kernel_type, degree
2152
2153            KernelType kernel_type = param.KernelType;
2154
2155            if (param.Degree < 0)
2156                return "degree of polynomial kernel < 0";
2157
2158            // cache_size,eps,C,nu,p,shrinking
2159
2160            if (param.CacheSize <= 0)
2161                return "cache_size <= 0";
2162
2163            if (param.EPS <= 0)
2164                return "eps <= 0";
2165
2166            if (param.Gamma == 0)
2167                param.Gamma = 1.0 / prob.MaxIndex;
2168
2169            if (svm_type == SvmType.C_SVC ||
2170               svm_type == SvmType.EPSILON_SVR ||
2171               svm_type == SvmType.NU_SVR)
2172                if (param.C <= 0)
2173                    return "C <= 0";
2174
2175            if (svm_type == SvmType.NU_SVC ||
2176               svm_type == SvmType.ONE_CLASS ||
2177               svm_type == SvmType.NU_SVR)
2178                if (param.Nu <= 0 || param.Nu > 1)
2179                    return "nu <= 0 or nu > 1";
2180
2181            if (svm_type == SvmType.EPSILON_SVR)
2182                if (param.P < 0)
2183                    return "p < 0";
2184
2185            if (param.Probability &&
2186               svm_type == SvmType.ONE_CLASS)
2187                return "one-class SVM probability output not supported yet";
2188
2189            // check whether nu-svc is feasible
2190
2191            if (svm_type == SvmType.NU_SVC)
2192            {
2193                int l = prob.Count;
2194                int Max_nr_class = 16;
2195                int nr_class = 0;
2196                int[] label = new int[Max_nr_class];
2197                int[] count = new int[Max_nr_class];
2198
2199                int i;
2200                for (i = 0; i < l; i++)
2201                {
2202                    int this_label = (int)prob.Y[i];
2203                    int j;
2204                    for (j = 0; j < nr_class; j++)
2205                        if (this_label == label[j])
2206                        {
2207                            ++count[j];
2208                            break;
2209                        }
2210
2211                    if (j == nr_class)
2212                    {
2213                        if (nr_class == Max_nr_class)
2214                        {
2215                            Max_nr_class *= 2;
2216                            int[] new_data = new int[Max_nr_class];
2217                            Array.Copy(label, 0, new_data, 0, label.Length);
2218                            label = new_data;
2219
2220                            new_data = new int[Max_nr_class];
2221                            Array.Copy(count, 0, new_data, 0, count.Length);
2222                            count = new_data;
2223                        }
2224                        label[nr_class] = this_label;
2225                        count[nr_class] = 1;
2226                        ++nr_class;
2227                    }
2228                }
2229
2230                for (i = 0; i < nr_class; i++)
2231                {
2232                    int n1 = count[i];
2233                    for (int j = i + 1; j < nr_class; j++)
2234                    {
2235                        int n2 = count[j];
2236                        if (param.Nu * (n1 + n2) / 2 > Math.Min(n1, n2))
2237                            return "specified nu is infeasible";
2238                    }
2239                }
2240            }
2241
2242            return null;
2243        }
2244
2245        public static int svm_check_probability_model(Model model)
2246        {
2247            if (((model.Parameter.SvmType == SvmType.C_SVC || model.Parameter.SvmType == SvmType.NU_SVC) &&
2248            model.PairwiseProbabilityA != null && model.PairwiseProbabilityB != null) ||
2249            ((model.Parameter.SvmType == SvmType.EPSILON_SVR || model.Parameter.SvmType == SvmType.NU_SVR) &&
2250             model.PairwiseProbabilityA != null))
2251                return 1;
2252            else
2253                return 0;
2254        }
2255    }
2256}
Note: See TracBrowser for help on using the repository browser.