Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/bdss.cs @ 2806

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

Added plugin for new version of ALGLIB. #875 (Update ALGLIB sources)

File size: 46.5 KB
Line 
1/*************************************************************************
2Copyright 2008 by Sergey Bochkanov (ALGLIB project).
3
4>>> SOURCE LICENSE >>>
5This program is free software; you can redistribute it and/or modify
6it under the terms of the GNU General Public License as published by
7the Free Software Foundation (www.fsf.org); either version 2 of the
8License, or (at your option) any later version.
9
10This program is distributed in the hope that it will be useful,
11but WITHOUT ANY WARRANTY; without even the implied warranty of
12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13GNU General Public License for more details.
14
15A copy of the GNU General Public License is available at
16http://www.fsf.org/licensing/licenses
17
18>>> END OF LICENSE >>>
19*************************************************************************/
20
21using System;
22
23namespace alglib
24{
25    public class bdss
26    {
27        public struct cvreport
28        {
29            public double relclserror;
30            public double avgce;
31            public double rmserror;
32            public double avgerror;
33            public double avgrelerror;
34        };
35
36
37
38
39        /*************************************************************************
40        This set of routines (DSErrAllocate, DSErrAccumulate, DSErrFinish)
41        calculates different error functions (classification error, cross-entropy,
42        rms, avg, avg.rel errors).
43
44        1. DSErrAllocate prepares buffer.
45        2. DSErrAccumulate accumulates individual errors:
46            * Y contains predicted output (posterior probabilities for classification)
47            * DesiredY contains desired output (class number for classification)
48        3. DSErrFinish outputs results:
49           * Buf[0] contains relative classification error (zero for regression tasks)
50           * Buf[1] contains avg. cross-entropy (zero for regression tasks)
51           * Buf[2] contains rms error (regression, classification)
52           * Buf[3] contains average error (regression, classification)
53           * Buf[4] contains average relative error (regression, classification)
54           
55        NOTES(1):
56            "NClasses>0" means that we have classification task.
57            "NClasses<0" means regression task with -NClasses real outputs.
58
59        NOTES(2):
60            rms. avg, avg.rel errors for classification tasks are interpreted as
61            errors in posterior probabilities with respect to probabilities given
62            by training/test set.
63
64          -- ALGLIB --
65             Copyright 11.01.2009 by Bochkanov Sergey
66        *************************************************************************/
67        public static void dserrallocate(int nclasses,
68            ref double[] buf)
69        {
70            buf = new double[7+1];
71            buf[0] = 0;
72            buf[1] = 0;
73            buf[2] = 0;
74            buf[3] = 0;
75            buf[4] = 0;
76            buf[5] = nclasses;
77            buf[6] = 0;
78            buf[7] = 0;
79        }
80
81
82        /*************************************************************************
83        See DSErrAllocate for comments on this routine.
84
85          -- ALGLIB --
86             Copyright 11.01.2009 by Bochkanov Sergey
87        *************************************************************************/
88        public static void dserraccumulate(ref double[] buf,
89            ref double[] y,
90            ref double[] desiredy)
91        {
92            int nclasses = 0;
93            int nout = 0;
94            int offs = 0;
95            int mmax = 0;
96            int rmax = 0;
97            int j = 0;
98            double v = 0;
99            double ev = 0;
100
101            offs = 5;
102            nclasses = (int)Math.Round(buf[offs]);
103            if( nclasses>0 )
104            {
105               
106                //
107                // Classification
108                //
109                rmax = (int)Math.Round(desiredy[0]);
110                mmax = 0;
111                for(j=1; j<=nclasses-1; j++)
112                {
113                    if( (double)(y[j])>(double)(y[mmax]) )
114                    {
115                        mmax = j;
116                    }
117                }
118                if( mmax!=rmax )
119                {
120                    buf[0] = buf[0]+1;
121                }
122                if( (double)(y[rmax])>(double)(0) )
123                {
124                    buf[1] = buf[1]-Math.Log(y[rmax]);
125                }
126                else
127                {
128                    buf[1] = buf[1]+Math.Log(AP.Math.MaxRealNumber);
129                }
130                for(j=0; j<=nclasses-1; j++)
131                {
132                    v = y[j];
133                    if( j==rmax )
134                    {
135                        ev = 1;
136                    }
137                    else
138                    {
139                        ev = 0;
140                    }
141                    buf[2] = buf[2]+AP.Math.Sqr(v-ev);
142                    buf[3] = buf[3]+Math.Abs(v-ev);
143                    if( (double)(ev)!=(double)(0) )
144                    {
145                        buf[4] = buf[4]+Math.Abs((v-ev)/ev);
146                        buf[offs+2] = buf[offs+2]+1;
147                    }
148                }
149                buf[offs+1] = buf[offs+1]+1;
150            }
151            else
152            {
153               
154                //
155                // Regression
156                //
157                nout = -nclasses;
158                rmax = 0;
159                for(j=1; j<=nout-1; j++)
160                {
161                    if( (double)(desiredy[j])>(double)(desiredy[rmax]) )
162                    {
163                        rmax = j;
164                    }
165                }
166                mmax = 0;
167                for(j=1; j<=nout-1; j++)
168                {
169                    if( (double)(y[j])>(double)(y[mmax]) )
170                    {
171                        mmax = j;
172                    }
173                }
174                if( mmax!=rmax )
175                {
176                    buf[0] = buf[0]+1;
177                }
178                for(j=0; j<=nout-1; j++)
179                {
180                    v = y[j];
181                    ev = desiredy[j];
182                    buf[2] = buf[2]+AP.Math.Sqr(v-ev);
183                    buf[3] = buf[3]+Math.Abs(v-ev);
184                    if( (double)(ev)!=(double)(0) )
185                    {
186                        buf[4] = buf[4]+Math.Abs((v-ev)/ev);
187                        buf[offs+2] = buf[offs+2]+1;
188                    }
189                }
190                buf[offs+1] = buf[offs+1]+1;
191            }
192        }
193
194
195        /*************************************************************************
196        See DSErrAllocate for comments on this routine.
197
198          -- ALGLIB --
199             Copyright 11.01.2009 by Bochkanov Sergey
200        *************************************************************************/
201        public static void dserrfinish(ref double[] buf)
202        {
203            int nout = 0;
204            int offs = 0;
205
206            offs = 5;
207            nout = Math.Abs((int)Math.Round(buf[offs]));
208            if( (double)(buf[offs+1])!=(double)(0) )
209            {
210                buf[0] = buf[0]/buf[offs+1];
211                buf[1] = buf[1]/buf[offs+1];
212                buf[2] = Math.Sqrt(buf[2]/(nout*buf[offs+1]));
213                buf[3] = buf[3]/(nout*buf[offs+1]);
214            }
215            if( (double)(buf[offs+2])!=(double)(0) )
216            {
217                buf[4] = buf[4]/buf[offs+2];
218            }
219        }
220
221
222        /*************************************************************************
223
224          -- ALGLIB --
225             Copyright 19.05.2008 by Bochkanov Sergey
226        *************************************************************************/
227        public static void dsnormalize(ref double[,] xy,
228            int npoints,
229            int nvars,
230            ref int info,
231            ref double[] means,
232            ref double[] sigmas)
233        {
234            int i = 0;
235            int j = 0;
236            double[] tmp = new double[0];
237            double mean = 0;
238            double variance = 0;
239            double skewness = 0;
240            double kurtosis = 0;
241            int i_ = 0;
242
243           
244            //
245            // Test parameters
246            //
247            if( npoints<=0 | nvars<1 )
248            {
249                info = -1;
250                return;
251            }
252            info = 1;
253           
254            //
255            // Standartization
256            //
257            means = new double[nvars-1+1];
258            sigmas = new double[nvars-1+1];
259            tmp = new double[npoints-1+1];
260            for(j=0; j<=nvars-1; j++)
261            {
262                for(i_=0; i_<=npoints-1;i_++)
263                {
264                    tmp[i_] = xy[i_,j];
265                }
266                descriptivestatistics.calculatemoments(ref tmp, npoints, ref mean, ref variance, ref skewness, ref kurtosis);
267                means[j] = mean;
268                sigmas[j] = Math.Sqrt(variance);
269                if( (double)(sigmas[j])==(double)(0) )
270                {
271                    sigmas[j] = 1;
272                }
273                for(i=0; i<=npoints-1; i++)
274                {
275                    xy[i,j] = (xy[i,j]-means[j])/sigmas[j];
276                }
277            }
278        }
279
280
281        /*************************************************************************
282
283          -- ALGLIB --
284             Copyright 19.05.2008 by Bochkanov Sergey
285        *************************************************************************/
286        public static void dsnormalizec(ref double[,] xy,
287            int npoints,
288            int nvars,
289            ref int info,
290            ref double[] means,
291            ref double[] sigmas)
292        {
293            int j = 0;
294            double[] tmp = new double[0];
295            double mean = 0;
296            double variance = 0;
297            double skewness = 0;
298            double kurtosis = 0;
299            int i_ = 0;
300
301           
302            //
303            // Test parameters
304            //
305            if( npoints<=0 | nvars<1 )
306            {
307                info = -1;
308                return;
309            }
310            info = 1;
311           
312            //
313            // Standartization
314            //
315            means = new double[nvars-1+1];
316            sigmas = new double[nvars-1+1];
317            tmp = new double[npoints-1+1];
318            for(j=0; j<=nvars-1; j++)
319            {
320                for(i_=0; i_<=npoints-1;i_++)
321                {
322                    tmp[i_] = xy[i_,j];
323                }
324                descriptivestatistics.calculatemoments(ref tmp, npoints, ref mean, ref variance, ref skewness, ref kurtosis);
325                means[j] = mean;
326                sigmas[j] = Math.Sqrt(variance);
327                if( (double)(sigmas[j])==(double)(0) )
328                {
329                    sigmas[j] = 1;
330                }
331            }
332        }
333
334
335        /*************************************************************************
336
337          -- ALGLIB --
338             Copyright 19.05.2008 by Bochkanov Sergey
339        *************************************************************************/
340        public static double dsgetmeanmindistance(ref double[,] xy,
341            int npoints,
342            int nvars)
343        {
344            double result = 0;
345            int i = 0;
346            int j = 0;
347            double[] tmp = new double[0];
348            double[] tmp2 = new double[0];
349            double v = 0;
350            int i_ = 0;
351
352           
353            //
354            // Test parameters
355            //
356            if( npoints<=0 | nvars<1 )
357            {
358                result = 0;
359                return result;
360            }
361           
362            //
363            // Process
364            //
365            tmp = new double[npoints-1+1];
366            for(i=0; i<=npoints-1; i++)
367            {
368                tmp[i] = AP.Math.MaxRealNumber;
369            }
370            tmp2 = new double[nvars-1+1];
371            for(i=0; i<=npoints-1; i++)
372            {
373                for(j=i+1; j<=npoints-1; j++)
374                {
375                    for(i_=0; i_<=nvars-1;i_++)
376                    {
377                        tmp2[i_] = xy[i,i_];
378                    }
379                    for(i_=0; i_<=nvars-1;i_++)
380                    {
381                        tmp2[i_] = tmp2[i_] - xy[j,i_];
382                    }
383                    v = 0.0;
384                    for(i_=0; i_<=nvars-1;i_++)
385                    {
386                        v += tmp2[i_]*tmp2[i_];
387                    }
388                    v = Math.Sqrt(v);
389                    tmp[i] = Math.Min(tmp[i], v);
390                    tmp[j] = Math.Min(tmp[j], v);
391                }
392            }
393            result = 0;
394            for(i=0; i<=npoints-1; i++)
395            {
396                result = result+tmp[i]/npoints;
397            }
398            return result;
399        }
400
401
402        /*************************************************************************
403
404          -- ALGLIB --
405             Copyright 19.05.2008 by Bochkanov Sergey
406        *************************************************************************/
407        public static void dstie(ref double[] a,
408            int n,
409            ref int[] ties,
410            ref int tiecount,
411            ref int[] p1,
412            ref int[] p2)
413        {
414            int i = 0;
415            int k = 0;
416            int[] tmp = new int[0];
417
418           
419            //
420            // Special case
421            //
422            if( n<=0 )
423            {
424                tiecount = 0;
425                return;
426            }
427           
428            //
429            // Sort A
430            //
431            tsort.tagsort(ref a, n, ref p1, ref p2);
432           
433            //
434            // Process ties
435            //
436            tiecount = 1;
437            for(i=1; i<=n-1; i++)
438            {
439                if( (double)(a[i])!=(double)(a[i-1]) )
440                {
441                    tiecount = tiecount+1;
442                }
443            }
444            ties = new int[tiecount+1];
445            ties[0] = 0;
446            k = 1;
447            for(i=1; i<=n-1; i++)
448            {
449                if( (double)(a[i])!=(double)(a[i-1]) )
450                {
451                    ties[k] = i;
452                    k = k+1;
453                }
454            }
455            ties[tiecount] = n;
456        }
457
458
459        /*************************************************************************
460
461          -- ALGLIB --
462             Copyright 11.12.2008 by Bochkanov Sergey
463        *************************************************************************/
464        public static void dstiefasti(ref double[] a,
465            ref int[] b,
466            int n,
467            ref int[] ties,
468            ref int tiecount)
469        {
470            int i = 0;
471            int k = 0;
472            int[] tmp = new int[0];
473
474           
475            //
476            // Special case
477            //
478            if( n<=0 )
479            {
480                tiecount = 0;
481                return;
482            }
483           
484            //
485            // Sort A
486            //
487            tsort.tagsortfasti(ref a, ref b, n);
488           
489            //
490            // Process ties
491            //
492            ties[0] = 0;
493            k = 1;
494            for(i=1; i<=n-1; i++)
495            {
496                if( (double)(a[i])!=(double)(a[i-1]) )
497                {
498                    ties[k] = i;
499                    k = k+1;
500                }
501            }
502            ties[k] = n;
503            tiecount = k;
504        }
505
506
507        /*************************************************************************
508        Optimal partition, internal subroutine.
509
510          -- ALGLIB --
511             Copyright 22.05.2008 by Bochkanov Sergey
512        *************************************************************************/
513        public static void dsoptimalsplit2(double[] a,
514            int[] c,
515            int n,
516            ref int info,
517            ref double threshold,
518            ref double pal,
519            ref double pbl,
520            ref double par,
521            ref double pbr,
522            ref double cve)
523        {
524            int i = 0;
525            int t = 0;
526            double s = 0;
527            int[] ties = new int[0];
528            int tiecount = 0;
529            int[] p1 = new int[0];
530            int[] p2 = new int[0];
531            int k = 0;
532            int koptimal = 0;
533            double pak = 0;
534            double pbk = 0;
535            double cvoptimal = 0;
536            double cv = 0;
537
538            a = (double[])a.Clone();
539            c = (int[])c.Clone();
540
541           
542            //
543            // Test for errors in inputs
544            //
545            if( n<=0 )
546            {
547                info = -1;
548                return;
549            }
550            for(i=0; i<=n-1; i++)
551            {
552                if( c[i]!=0 & c[i]!=1 )
553                {
554                    info = -2;
555                    return;
556                }
557            }
558            info = 1;
559           
560            //
561            // Tie
562            //
563            dstie(ref a, n, ref ties, ref tiecount, ref p1, ref p2);
564            for(i=0; i<=n-1; i++)
565            {
566                if( p2[i]!=i )
567                {
568                    t = c[i];
569                    c[i] = c[p2[i]];
570                    c[p2[i]] = t;
571                }
572            }
573           
574            //
575            // Special case: number of ties is 1.
576            //
577            // NOTE: we assume that P[i,j] equals to 0 or 1,
578            //       intermediate values are not allowed.
579            //
580            if( tiecount==1 )
581            {
582                info = -3;
583                return;
584            }
585           
586            //
587            // General case, number of ties > 1
588            //
589            // NOTE: we assume that P[i,j] equals to 0 or 1,
590            //       intermediate values are not allowed.
591            //
592            pal = 0;
593            pbl = 0;
594            par = 0;
595            pbr = 0;
596            for(i=0; i<=n-1; i++)
597            {
598                if( c[i]==0 )
599                {
600                    par = par+1;
601                }
602                if( c[i]==1 )
603                {
604                    pbr = pbr+1;
605                }
606            }
607            koptimal = -1;
608            cvoptimal = AP.Math.MaxRealNumber;
609            for(k=0; k<=tiecount-2; k++)
610            {
611               
612                //
613                // first, obtain information about K-th tie which is
614                // moved from R-part to L-part
615                //
616                pak = 0;
617                pbk = 0;
618                for(i=ties[k]; i<=ties[k+1]-1; i++)
619                {
620                    if( c[i]==0 )
621                    {
622                        pak = pak+1;
623                    }
624                    if( c[i]==1 )
625                    {
626                        pbk = pbk+1;
627                    }
628                }
629               
630                //
631                // Calculate cross-validation CE
632                //
633                cv = 0;
634                cv = cv-xlny(pal+pak, (pal+pak)/(pal+pak+pbl+pbk+1));
635                cv = cv-xlny(pbl+pbk, (pbl+pbk)/(pal+pak+1+pbl+pbk));
636                cv = cv-xlny(par-pak, (par-pak)/(par-pak+pbr-pbk+1));
637                cv = cv-xlny(pbr-pbk, (pbr-pbk)/(par-pak+1+pbr-pbk));
638               
639                //
640                // Compare with best
641                //
642                if( (double)(cv)<(double)(cvoptimal) )
643                {
644                    cvoptimal = cv;
645                    koptimal = k;
646                }
647               
648                //
649                // update
650                //
651                pal = pal+pak;
652                pbl = pbl+pbk;
653                par = par-pak;
654                pbr = pbr-pbk;
655            }
656            cve = cvoptimal;
657            threshold = 0.5*(a[ties[koptimal]]+a[ties[koptimal+1]]);
658            pal = 0;
659            pbl = 0;
660            par = 0;
661            pbr = 0;
662            for(i=0; i<=n-1; i++)
663            {
664                if( (double)(a[i])<(double)(threshold) )
665                {
666                    if( c[i]==0 )
667                    {
668                        pal = pal+1;
669                    }
670                    else
671                    {
672                        pbl = pbl+1;
673                    }
674                }
675                else
676                {
677                    if( c[i]==0 )
678                    {
679                        par = par+1;
680                    }
681                    else
682                    {
683                        pbr = pbr+1;
684                    }
685                }
686            }
687            s = pal+pbl;
688            pal = pal/s;
689            pbl = pbl/s;
690            s = par+pbr;
691            par = par/s;
692            pbr = pbr/s;
693        }
694
695
696        /*************************************************************************
697        Optimal partition, internal subroutine. Fast version.
698
699        Accepts:
700            A       array[0..N-1]       array of attributes     array[0..N-1]
701            C       array[0..N-1]       array of class labels
702            TiesBuf array[0..N]         temporaries (ties)
703            CntBuf  array[0..2*NC-1]    temporaries (counts)
704            Alpha                       centering factor (0<=alpha<=1, recommended value - 0.05)
705           
706        Output:
707            Info    error code (">0"=OK, "<0"=bad)
708            RMS     training set RMS error
709            CVRMS   leave-one-out RMS error
710           
711        Note:
712            content of all arrays is changed by subroutine
713
714          -- ALGLIB --
715             Copyright 11.12.2008 by Bochkanov Sergey
716        *************************************************************************/
717        public static void dsoptimalsplit2fast(ref double[] a,
718            ref int[] c,
719            ref int[] tiesbuf,
720            ref int[] cntbuf,
721            int n,
722            int nc,
723            double alpha,
724            ref int info,
725            ref double threshold,
726            ref double rms,
727            ref double cvrms)
728        {
729            int i = 0;
730            int k = 0;
731            int cl = 0;
732            int tiecount = 0;
733            double cbest = 0;
734            double cc = 0;
735            int koptimal = 0;
736            int sl = 0;
737            int sr = 0;
738            double v = 0;
739            double w = 0;
740            double x = 0;
741
742           
743            //
744            // Test for errors in inputs
745            //
746            if( n<=0 | nc<2 )
747            {
748                info = -1;
749                return;
750            }
751            for(i=0; i<=n-1; i++)
752            {
753                if( c[i]<0 | c[i]>=nc )
754                {
755                    info = -2;
756                    return;
757                }
758            }
759            info = 1;
760           
761            //
762            // Tie
763            //
764            dstiefasti(ref a, ref c, n, ref tiesbuf, ref tiecount);
765           
766            //
767            // Special case: number of ties is 1.
768            //
769            if( tiecount==1 )
770            {
771                info = -3;
772                return;
773            }
774           
775            //
776            // General case, number of ties > 1
777            //
778            for(i=0; i<=2*nc-1; i++)
779            {
780                cntbuf[i] = 0;
781            }
782            for(i=0; i<=n-1; i++)
783            {
784                cntbuf[nc+c[i]] = cntbuf[nc+c[i]]+1;
785            }
786            koptimal = -1;
787            threshold = a[n-1];
788            cbest = AP.Math.MaxRealNumber;
789            sl = 0;
790            sr = n;
791            for(k=0; k<=tiecount-2; k++)
792            {
793               
794                //
795                // first, move Kth tie from right to left
796                //
797                for(i=tiesbuf[k]; i<=tiesbuf[k+1]-1; i++)
798                {
799                    cl = c[i];
800                    cntbuf[cl] = cntbuf[cl]+1;
801                    cntbuf[nc+cl] = cntbuf[nc+cl]-1;
802                }
803                sl = sl+(tiesbuf[k+1]-tiesbuf[k]);
804                sr = sr-(tiesbuf[k+1]-tiesbuf[k]);
805               
806                //
807                // Calculate RMS error
808                //
809                v = 0;
810                for(i=0; i<=nc-1; i++)
811                {
812                    w = cntbuf[i];
813                    v = v+w*AP.Math.Sqr(w/sl-1);
814                    v = v+(sl-w)*AP.Math.Sqr(w/sl);
815                    w = cntbuf[nc+i];
816                    v = v+w*AP.Math.Sqr(w/sr-1);
817                    v = v+(sr-w)*AP.Math.Sqr(w/sr);
818                }
819                v = Math.Sqrt(v/(nc*n));
820               
821                //
822                // Compare with best
823                //
824                x = (double)(2*sl)/((double)(sl+sr))-1;
825                cc = v*(1-alpha+alpha*AP.Math.Sqr(x));
826                if( (double)(cc)<(double)(cbest) )
827                {
828                   
829                    //
830                    // store split
831                    //
832                    rms = v;
833                    koptimal = k;
834                    cbest = cc;
835                   
836                    //
837                    // calculate CVRMS error
838                    //
839                    cvrms = 0;
840                    for(i=0; i<=nc-1; i++)
841                    {
842                        if( sl>1 )
843                        {
844                            w = cntbuf[i];
845                            cvrms = cvrms+w*AP.Math.Sqr((w-1)/(sl-1)-1);
846                            cvrms = cvrms+(sl-w)*AP.Math.Sqr(w/(sl-1));
847                        }
848                        else
849                        {
850                            w = cntbuf[i];
851                            cvrms = cvrms+w*AP.Math.Sqr((double)(1)/(double)(nc)-1);
852                            cvrms = cvrms+(sl-w)*AP.Math.Sqr((double)(1)/(double)(nc));
853                        }
854                        if( sr>1 )
855                        {
856                            w = cntbuf[nc+i];
857                            cvrms = cvrms+w*AP.Math.Sqr((w-1)/(sr-1)-1);
858                            cvrms = cvrms+(sr-w)*AP.Math.Sqr(w/(sr-1));
859                        }
860                        else
861                        {
862                            w = cntbuf[nc+i];
863                            cvrms = cvrms+w*AP.Math.Sqr((double)(1)/(double)(nc)-1);
864                            cvrms = cvrms+(sr-w)*AP.Math.Sqr((double)(1)/(double)(nc));
865                        }
866                    }
867                    cvrms = Math.Sqrt(cvrms/(nc*n));
868                }
869            }
870           
871            //
872            // Calculate threshold.
873            // Code is a bit complicated because there can be such
874            // numbers that 0.5(A+B) equals to A or B (if A-B=epsilon)
875            //
876            threshold = 0.5*(a[tiesbuf[koptimal]]+a[tiesbuf[koptimal+1]]);
877            if( (double)(threshold)<=(double)(a[tiesbuf[koptimal]]) )
878            {
879                threshold = a[tiesbuf[koptimal+1]];
880            }
881        }
882
883
884        /*************************************************************************
885        Automatic non-optimal discretization, internal subroutine.
886
887          -- ALGLIB --
888             Copyright 22.05.2008 by Bochkanov Sergey
889        *************************************************************************/
890        public static void dssplitk(double[] a,
891            int[] c,
892            int n,
893            int nc,
894            int kmax,
895            ref int info,
896            ref double[] thresholds,
897            ref int ni,
898            ref double cve)
899        {
900            int i = 0;
901            int j = 0;
902            int j1 = 0;
903            int k = 0;
904            int[] ties = new int[0];
905            int tiecount = 0;
906            int[] p1 = new int[0];
907            int[] p2 = new int[0];
908            int[] cnt = new int[0];
909            double v2 = 0;
910            int bestk = 0;
911            double bestcve = 0;
912            int[] bestsizes = new int[0];
913            double curcve = 0;
914            int[] cursizes = new int[0];
915
916            a = (double[])a.Clone();
917            c = (int[])c.Clone();
918
919           
920            //
921            // Test for errors in inputs
922            //
923            if( n<=0 | nc<2 | kmax<2 )
924            {
925                info = -1;
926                return;
927            }
928            for(i=0; i<=n-1; i++)
929            {
930                if( c[i]<0 | c[i]>=nc )
931                {
932                    info = -2;
933                    return;
934                }
935            }
936            info = 1;
937           
938            //
939            // Tie
940            //
941            dstie(ref a, n, ref ties, ref tiecount, ref p1, ref p2);
942            for(i=0; i<=n-1; i++)
943            {
944                if( p2[i]!=i )
945                {
946                    k = c[i];
947                    c[i] = c[p2[i]];
948                    c[p2[i]] = k;
949                }
950            }
951           
952            //
953            // Special cases
954            //
955            if( tiecount==1 )
956            {
957                info = -3;
958                return;
959            }
960           
961            //
962            // General case:
963            // 0. allocate arrays
964            //
965            kmax = Math.Min(kmax, tiecount);
966            bestsizes = new int[kmax-1+1];
967            cursizes = new int[kmax-1+1];
968            cnt = new int[nc-1+1];
969           
970            //
971            // General case:
972            // 1. prepare "weak" solution (two subintervals, divided at median)
973            //
974            v2 = AP.Math.MaxRealNumber;
975            j = -1;
976            for(i=1; i<=tiecount-1; i++)
977            {
978                if( (double)(Math.Abs(ties[i]-0.5*(n-1)))<(double)(v2) )
979                {
980                    v2 = Math.Abs(ties[i]-0.5*n);
981                    j = i;
982                }
983            }
984            System.Diagnostics.Debug.Assert(j>0, "DSSplitK: internal error #1!");
985            bestk = 2;
986            bestsizes[0] = ties[j];
987            bestsizes[1] = n-j;
988            bestcve = 0;
989            for(i=0; i<=nc-1; i++)
990            {
991                cnt[i] = 0;
992            }
993            for(i=0; i<=j-1; i++)
994            {
995                tieaddc(ref c, ref ties, i, nc, ref cnt);
996            }
997            bestcve = bestcve+getcv(ref cnt, nc);
998            for(i=0; i<=nc-1; i++)
999            {
1000                cnt[i] = 0;
1001            }
1002            for(i=j; i<=tiecount-1; i++)
1003            {
1004                tieaddc(ref c, ref ties, i, nc, ref cnt);
1005            }
1006            bestcve = bestcve+getcv(ref cnt, nc);
1007           
1008            //
1009            // General case:
1010            // 2. Use greedy algorithm to find sub-optimal split in O(KMax*N) time
1011            //
1012            for(k=2; k<=kmax; k++)
1013            {
1014               
1015                //
1016                // Prepare greedy K-interval split
1017                //
1018                for(i=0; i<=k-1; i++)
1019                {
1020                    cursizes[i] = 0;
1021                }
1022                i = 0;
1023                j = 0;
1024                while( j<=tiecount-1 & i<=k-1 )
1025                {
1026                   
1027                    //
1028                    // Rule: I-th bin is empty, fill it
1029                    //
1030                    if( cursizes[i]==0 )
1031                    {
1032                        cursizes[i] = ties[j+1]-ties[j];
1033                        j = j+1;
1034                        continue;
1035                    }
1036                   
1037                    //
1038                    // Rule: (K-1-I) bins left, (K-1-I) ties left (1 tie per bin); next bin
1039                    //
1040                    if( tiecount-j==k-1-i )
1041                    {
1042                        i = i+1;
1043                        continue;
1044                    }
1045                   
1046                    //
1047                    // Rule: last bin, always place in current
1048                    //
1049                    if( i==k-1 )
1050                    {
1051                        cursizes[i] = cursizes[i]+ties[j+1]-ties[j];
1052                        j = j+1;
1053                        continue;
1054                    }
1055                   
1056                    //
1057                    // Place J-th tie in I-th bin, or leave for I+1-th bin.
1058                    //
1059                    if( (double)(Math.Abs(cursizes[i]+ties[j+1]-ties[j]-(double)(n)/(double)(k)))<(double)(Math.Abs(cursizes[i]-(double)(n)/(double)(k))) )
1060                    {
1061                        cursizes[i] = cursizes[i]+ties[j+1]-ties[j];
1062                        j = j+1;
1063                    }
1064                    else
1065                    {
1066                        i = i+1;
1067                    }
1068                }
1069                System.Diagnostics.Debug.Assert(cursizes[k-1]!=0 & j==tiecount, "DSSplitK: internal error #1");
1070               
1071                //
1072                // Calculate CVE
1073                //
1074                curcve = 0;
1075                j = 0;
1076                for(i=0; i<=k-1; i++)
1077                {
1078                    for(j1=0; j1<=nc-1; j1++)
1079                    {
1080                        cnt[j1] = 0;
1081                    }
1082                    for(j1=j; j1<=j+cursizes[i]-1; j1++)
1083                    {
1084                        cnt[c[j1]] = cnt[c[j1]]+1;
1085                    }
1086                    curcve = curcve+getcv(ref cnt, nc);
1087                    j = j+cursizes[i];
1088                }
1089               
1090                //
1091                // Choose best variant
1092                //
1093                if( (double)(curcve)<(double)(bestcve) )
1094                {
1095                    for(i=0; i<=k-1; i++)
1096                    {
1097                        bestsizes[i] = cursizes[i];
1098                    }
1099                    bestcve = curcve;
1100                    bestk = k;
1101                }
1102            }
1103           
1104            //
1105            // Transform from sizes to thresholds
1106            //
1107            cve = bestcve;
1108            ni = bestk;
1109            thresholds = new double[ni-2+1];
1110            j = bestsizes[0];
1111            for(i=1; i<=bestk-1; i++)
1112            {
1113                thresholds[i-1] = 0.5*(a[j-1]+a[j]);
1114                j = j+bestsizes[i];
1115            }
1116        }
1117
1118
1119        /*************************************************************************
1120        Automatic optimal discretization, internal subroutine.
1121
1122          -- ALGLIB --
1123             Copyright 22.05.2008 by Bochkanov Sergey
1124        *************************************************************************/
1125        public static void dsoptimalsplitk(double[] a,
1126            int[] c,
1127            int n,
1128            int nc,
1129            int kmax,
1130            ref int info,
1131            ref double[] thresholds,
1132            ref int ni,
1133            ref double cve)
1134        {
1135            int i = 0;
1136            int j = 0;
1137            int s = 0;
1138            int jl = 0;
1139            int jr = 0;
1140            double v2 = 0;
1141            int[] ties = new int[0];
1142            int tiecount = 0;
1143            int[] p1 = new int[0];
1144            int[] p2 = new int[0];
1145            double cvtemp = 0;
1146            int[] cnt = new int[0];
1147            int[] cnt2 = new int[0];
1148            double[,] cv = new double[0,0];
1149            int[,] splits = new int[0,0];
1150            int k = 0;
1151            int koptimal = 0;
1152            double cvoptimal = 0;
1153
1154            a = (double[])a.Clone();
1155            c = (int[])c.Clone();
1156
1157           
1158            //
1159            // Test for errors in inputs
1160            //
1161            if( n<=0 | nc<2 | kmax<2 )
1162            {
1163                info = -1;
1164                return;
1165            }
1166            for(i=0; i<=n-1; i++)
1167            {
1168                if( c[i]<0 | c[i]>=nc )
1169                {
1170                    info = -2;
1171                    return;
1172                }
1173            }
1174            info = 1;
1175           
1176            //
1177            // Tie
1178            //
1179            dstie(ref a, n, ref ties, ref tiecount, ref p1, ref p2);
1180            for(i=0; i<=n-1; i++)
1181            {
1182                if( p2[i]!=i )
1183                {
1184                    k = c[i];
1185                    c[i] = c[p2[i]];
1186                    c[p2[i]] = k;
1187                }
1188            }
1189           
1190            //
1191            // Special cases
1192            //
1193            if( tiecount==1 )
1194            {
1195                info = -3;
1196                return;
1197            }
1198           
1199            //
1200            // General case
1201            // Use dynamic programming to find best split in O(KMax*NC*TieCount^2) time
1202            //
1203            kmax = Math.Min(kmax, tiecount);
1204            cv = new double[kmax-1+1, tiecount-1+1];
1205            splits = new int[kmax-1+1, tiecount-1+1];
1206            cnt = new int[nc-1+1];
1207            cnt2 = new int[nc-1+1];
1208            for(j=0; j<=nc-1; j++)
1209            {
1210                cnt[j] = 0;
1211            }
1212            for(j=0; j<=tiecount-1; j++)
1213            {
1214                tieaddc(ref c, ref ties, j, nc, ref cnt);
1215                splits[0,j] = 0;
1216                cv[0,j] = getcv(ref cnt, nc);
1217            }
1218            for(k=1; k<=kmax-1; k++)
1219            {
1220                for(j=0; j<=nc-1; j++)
1221                {
1222                    cnt[j] = 0;
1223                }
1224               
1225                //
1226                // Subtask size J in [K..TieCount-1]:
1227                // optimal K-splitting on ties from 0-th to J-th.
1228                //
1229                for(j=k; j<=tiecount-1; j++)
1230                {
1231                   
1232                    //
1233                    // Update Cnt - let it contain classes of ties from K-th to J-th
1234                    //
1235                    tieaddc(ref c, ref ties, j, nc, ref cnt);
1236                   
1237                    //
1238                    // Search for optimal split point S in [K..J]
1239                    //
1240                    for(i=0; i<=nc-1; i++)
1241                    {
1242                        cnt2[i] = cnt[i];
1243                    }
1244                    cv[k,j] = cv[k-1,j-1]+getcv(ref cnt2, nc);
1245                    splits[k,j] = j;
1246                    for(s=k+1; s<=j; s++)
1247                    {
1248                       
1249                        //
1250                        // Update Cnt2 - let it contain classes of ties from S-th to J-th
1251                        //
1252                        tiesubc(ref c, ref ties, s-1, nc, ref cnt2);
1253                       
1254                        //
1255                        // Calculate CVE
1256                        //
1257                        cvtemp = cv[k-1,s-1]+getcv(ref cnt2, nc);
1258                        if( (double)(cvtemp)<(double)(cv[k,j]) )
1259                        {
1260                            cv[k,j] = cvtemp;
1261                            splits[k,j] = s;
1262                        }
1263                    }
1264                }
1265            }
1266           
1267            //
1268            // Choose best partition, output result
1269            //
1270            koptimal = -1;
1271            cvoptimal = AP.Math.MaxRealNumber;
1272            for(k=0; k<=kmax-1; k++)
1273            {
1274                if( (double)(cv[k,tiecount-1])<(double)(cvoptimal) )
1275                {
1276                    cvoptimal = cv[k,tiecount-1];
1277                    koptimal = k;
1278                }
1279            }
1280            System.Diagnostics.Debug.Assert(koptimal>=0, "DSOptimalSplitK: internal error #1!");
1281            if( koptimal==0 )
1282            {
1283               
1284                //
1285                // Special case: best partition is one big interval.
1286                // Even 2-partition is not better.
1287                // This is possible when dealing with "weak" predictor variables.
1288                //
1289                // Make binary split as close to the median as possible.
1290                //
1291                v2 = AP.Math.MaxRealNumber;
1292                j = -1;
1293                for(i=1; i<=tiecount-1; i++)
1294                {
1295                    if( (double)(Math.Abs(ties[i]-0.5*(n-1)))<(double)(v2) )
1296                    {
1297                        v2 = Math.Abs(ties[i]-0.5*(n-1));
1298                        j = i;
1299                    }
1300                }
1301                System.Diagnostics.Debug.Assert(j>0, "DSOptimalSplitK: internal error #2!");
1302                thresholds = new double[0+1];
1303                thresholds[0] = 0.5*(a[ties[j-1]]+a[ties[j]]);
1304                ni = 2;
1305                cve = 0;
1306                for(i=0; i<=nc-1; i++)
1307                {
1308                    cnt[i] = 0;
1309                }
1310                for(i=0; i<=j-1; i++)
1311                {
1312                    tieaddc(ref c, ref ties, i, nc, ref cnt);
1313                }
1314                cve = cve+getcv(ref cnt, nc);
1315                for(i=0; i<=nc-1; i++)
1316                {
1317                    cnt[i] = 0;
1318                }
1319                for(i=j; i<=tiecount-1; i++)
1320                {
1321                    tieaddc(ref c, ref ties, i, nc, ref cnt);
1322                }
1323                cve = cve+getcv(ref cnt, nc);
1324            }
1325            else
1326            {
1327               
1328                //
1329                // General case: 2 or more intervals
1330                //
1331                thresholds = new double[koptimal-1+1];
1332                ni = koptimal+1;
1333                cve = cv[koptimal,tiecount-1];
1334                jl = splits[koptimal,tiecount-1];
1335                jr = tiecount-1;
1336                for(k=koptimal; k>=1; k--)
1337                {
1338                    thresholds[k-1] = 0.5*(a[ties[jl-1]]+a[ties[jl]]);
1339                    jr = jl-1;
1340                    jl = splits[k-1,jl-1];
1341                }
1342            }
1343        }
1344
1345
1346        /*************************************************************************
1347        Subroutine prepares K-fold split of the training set.
1348
1349        NOTES:
1350            "NClasses>0" means that we have classification task.
1351            "NClasses<0" means regression task with -NClasses real outputs.
1352
1353          -- ALGLIB --
1354             Copyright 11.01.2009 by Bochkanov Sergey
1355        *************************************************************************/
1356        private static void dskfoldsplit(ref double[,] xy,
1357            int npoints,
1358            int nclasses,
1359            int foldscount,
1360            bool stratifiedsplits,
1361            ref int[] folds)
1362        {
1363            int i = 0;
1364            int j = 0;
1365            int k = 0;
1366
1367           
1368            //
1369            // test parameters
1370            //
1371            System.Diagnostics.Debug.Assert(npoints>0, "DSKFoldSplit: wrong NPoints!");
1372            System.Diagnostics.Debug.Assert(nclasses>1 | nclasses<0, "DSKFoldSplit: wrong NClasses!");
1373            System.Diagnostics.Debug.Assert(foldscount>=2 & foldscount<=npoints, "DSKFoldSplit: wrong FoldsCount!");
1374            System.Diagnostics.Debug.Assert(!stratifiedsplits, "DSKFoldSplit: stratified splits are not supported!");
1375           
1376            //
1377            // Folds
1378            //
1379            folds = new int[npoints-1+1];
1380            for(i=0; i<=npoints-1; i++)
1381            {
1382                folds[i] = i*foldscount/npoints;
1383            }
1384            for(i=0; i<=npoints-2; i++)
1385            {
1386                j = i+AP.Math.RandomInteger(npoints-i);
1387                if( j!=i )
1388                {
1389                    k = folds[i];
1390                    folds[i] = folds[j];
1391                    folds[j] = k;
1392                }
1393            }
1394        }
1395
1396
1397        /*************************************************************************
1398        Internal function
1399        *************************************************************************/
1400        private static double xlny(double x,
1401            double y)
1402        {
1403            double result = 0;
1404
1405            if( (double)(x)==(double)(0) )
1406            {
1407                result = 0;
1408            }
1409            else
1410            {
1411                result = x*Math.Log(y);
1412            }
1413            return result;
1414        }
1415
1416
1417        /*************************************************************************
1418        Internal function,
1419        returns number of samples of class I in Cnt[I]
1420        *************************************************************************/
1421        private static double getcv(ref int[] cnt,
1422            int nc)
1423        {
1424            double result = 0;
1425            int i = 0;
1426            double s = 0;
1427
1428            s = 0;
1429            for(i=0; i<=nc-1; i++)
1430            {
1431                s = s+cnt[i];
1432            }
1433            result = 0;
1434            for(i=0; i<=nc-1; i++)
1435            {
1436                result = result-xlny(cnt[i], cnt[i]/(s+nc-1));
1437            }
1438            return result;
1439        }
1440
1441
1442        /*************************************************************************
1443        Internal function, adds number of samples of class I in tie NTie to Cnt[I]
1444        *************************************************************************/
1445        private static void tieaddc(ref int[] c,
1446            ref int[] ties,
1447            int ntie,
1448            int nc,
1449            ref int[] cnt)
1450        {
1451            int i = 0;
1452
1453            for(i=ties[ntie]; i<=ties[ntie+1]-1; i++)
1454            {
1455                cnt[c[i]] = cnt[c[i]]+1;
1456            }
1457        }
1458
1459
1460        /*************************************************************************
1461        Internal function, subtracts number of samples of class I in tie NTie to Cnt[I]
1462        *************************************************************************/
1463        private static void tiesubc(ref int[] c,
1464            ref int[] ties,
1465            int ntie,
1466            int nc,
1467            ref int[] cnt)
1468        {
1469            int i = 0;
1470
1471            for(i=ties[ntie]; i<=ties[ntie+1]-1; i++)
1472            {
1473                cnt[c[i]] = cnt[c[i]]-1;
1474            }
1475        }
1476
1477
1478        /*************************************************************************
1479        Internal function,
1480        returns number of samples of class I in Cnt[I]
1481        *************************************************************************/
1482        private static void tiegetc(ref int[] c,
1483            ref int[] ties,
1484            int ntie,
1485            int nc,
1486            ref int[] cnt)
1487        {
1488            int i = 0;
1489
1490            for(i=0; i<=nc-1; i++)
1491            {
1492                cnt[i] = 0;
1493            }
1494            for(i=ties[ntie]; i<=ties[ntie+1]-1; i++)
1495            {
1496                cnt[c[i]] = cnt[c[i]]+1;
1497            }
1498        }
1499    }
1500}
Note: See TracBrowser for help on using the repository browser.