Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/bdss.cs @ 2636

Last change on this file since 2636 was 2563, checked in by gkronber, 15 years ago

Updated ALGLIB to latest version. #751 (Plugin for for data-modeling with ANN (integrated into CEDMA))

File size: 46.7 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 i = 0;
294            int j = 0;
295            double[] tmp = new double[0];
296            double mean = 0;
297            double variance = 0;
298            double skewness = 0;
299            double kurtosis = 0;
300            int i_ = 0;
301
302           
303            //
304            // Test parameters
305            //
306            if( npoints<=0 | nvars<1 )
307            {
308                info = -1;
309                return;
310            }
311            info = 1;
312           
313            //
314            // Standartization
315            //
316            means = new double[nvars-1+1];
317            sigmas = new double[nvars-1+1];
318            tmp = new double[npoints-1+1];
319            for(j=0; j<=nvars-1; j++)
320            {
321                for(i_=0; i_<=npoints-1;i_++)
322                {
323                    tmp[i_] = xy[i_,j];
324                }
325                descriptivestatistics.calculatemoments(ref tmp, npoints, ref mean, ref variance, ref skewness, ref kurtosis);
326                means[j] = mean;
327                sigmas[j] = Math.Sqrt(variance);
328                if( (double)(sigmas[j])==(double)(0) )
329                {
330                    sigmas[j] = 1;
331                }
332            }
333        }
334
335
336        /*************************************************************************
337
338          -- ALGLIB --
339             Copyright 19.05.2008 by Bochkanov Sergey
340        *************************************************************************/
341        public static double dsgetmeanmindistance(ref double[,] xy,
342            int npoints,
343            int nvars)
344        {
345            double result = 0;
346            int i = 0;
347            int j = 0;
348            double[] tmp = new double[0];
349            double[] tmp2 = new double[0];
350            double v = 0;
351            int i_ = 0;
352
353           
354            //
355            // Test parameters
356            //
357            if( npoints<=0 | nvars<1 )
358            {
359                result = 0;
360                return result;
361            }
362           
363            //
364            // Process
365            //
366            tmp = new double[npoints-1+1];
367            for(i=0; i<=npoints-1; i++)
368            {
369                tmp[i] = AP.Math.MaxRealNumber;
370            }
371            tmp2 = new double[nvars-1+1];
372            for(i=0; i<=npoints-1; i++)
373            {
374                for(j=i+1; j<=npoints-1; j++)
375                {
376                    for(i_=0; i_<=nvars-1;i_++)
377                    {
378                        tmp2[i_] = xy[i,i_];
379                    }
380                    for(i_=0; i_<=nvars-1;i_++)
381                    {
382                        tmp2[i_] = tmp2[i_] - xy[j,i_];
383                    }
384                    v = 0.0;
385                    for(i_=0; i_<=nvars-1;i_++)
386                    {
387                        v += tmp2[i_]*tmp2[i_];
388                    }
389                    v = Math.Sqrt(v);
390                    tmp[i] = Math.Min(tmp[i], v);
391                    tmp[j] = Math.Min(tmp[j], v);
392                }
393            }
394            result = 0;
395            for(i=0; i<=npoints-1; i++)
396            {
397                result = result+tmp[i]/npoints;
398            }
399            return result;
400        }
401
402
403        /*************************************************************************
404
405          -- ALGLIB --
406             Copyright 19.05.2008 by Bochkanov Sergey
407        *************************************************************************/
408        public static void dstie(ref double[] a,
409            int n,
410            ref int[] ties,
411            ref int tiecount,
412            ref int[] p1,
413            ref int[] p2)
414        {
415            int i = 0;
416            int k = 0;
417            int[] tmp = new int[0];
418
419           
420            //
421            // Special case
422            //
423            if( n<=0 )
424            {
425                tiecount = 0;
426                return;
427            }
428           
429            //
430            // Sort A
431            //
432            tsort.tagsort(ref a, n, ref p1, ref p2);
433           
434            //
435            // Process ties
436            //
437            tiecount = 1;
438            for(i=1; i<=n-1; i++)
439            {
440                if( (double)(a[i])!=(double)(a[i-1]) )
441                {
442                    tiecount = tiecount+1;
443                }
444            }
445            ties = new int[tiecount+1];
446            ties[0] = 0;
447            k = 1;
448            for(i=1; i<=n-1; i++)
449            {
450                if( (double)(a[i])!=(double)(a[i-1]) )
451                {
452                    ties[k] = i;
453                    k = k+1;
454                }
455            }
456            ties[tiecount] = n;
457        }
458
459
460        /*************************************************************************
461
462          -- ALGLIB --
463             Copyright 11.12.2008 by Bochkanov Sergey
464        *************************************************************************/
465        public static void dstiefasti(ref double[] a,
466            ref int[] b,
467            int n,
468            ref int[] ties,
469            ref int tiecount)
470        {
471            int i = 0;
472            int k = 0;
473            int[] tmp = new int[0];
474
475           
476            //
477            // Special case
478            //
479            if( n<=0 )
480            {
481                tiecount = 0;
482                return;
483            }
484           
485            //
486            // Sort A
487            //
488            tsort.tagsortfasti(ref a, ref b, n);
489           
490            //
491            // Process ties
492            //
493            ties[0] = 0;
494            k = 1;
495            for(i=1; i<=n-1; i++)
496            {
497                if( (double)(a[i])!=(double)(a[i-1]) )
498                {
499                    ties[k] = i;
500                    k = k+1;
501                }
502            }
503            ties[k] = n;
504            tiecount = k;
505        }
506
507
508        /*************************************************************************
509        Optimal partition, internal subroutine.
510
511          -- ALGLIB --
512             Copyright 22.05.2008 by Bochkanov Sergey
513        *************************************************************************/
514        public static void dsoptimalsplit2(double[] a,
515            int[] c,
516            int n,
517            ref int info,
518            ref double threshold,
519            ref double pal,
520            ref double pbl,
521            ref double par,
522            ref double pbr,
523            ref double cve)
524        {
525            int i = 0;
526            int t = 0;
527            double s = 0;
528            double pea = 0;
529            double peb = 0;
530            int[] ties = new int[0];
531            int tiecount = 0;
532            int[] p1 = new int[0];
533            int[] p2 = new int[0];
534            double v1 = 0;
535            double v2 = 0;
536            int k = 0;
537            int koptimal = 0;
538            double pak = 0;
539            double pbk = 0;
540            double cvoptimal = 0;
541            double cv = 0;
542
543            a = (double[])a.Clone();
544            c = (int[])c.Clone();
545
546           
547            //
548            // Test for errors in inputs
549            //
550            if( n<=0 )
551            {
552                info = -1;
553                return;
554            }
555            for(i=0; i<=n-1; i++)
556            {
557                if( c[i]!=0 & c[i]!=1 )
558                {
559                    info = -2;
560                    return;
561                }
562            }
563            info = 1;
564           
565            //
566            // Tie
567            //
568            dstie(ref a, n, ref ties, ref tiecount, ref p1, ref p2);
569            for(i=0; i<=n-1; i++)
570            {
571                if( p2[i]!=i )
572                {
573                    t = c[i];
574                    c[i] = c[p2[i]];
575                    c[p2[i]] = t;
576                }
577            }
578           
579            //
580            // Special case: number of ties is 1.
581            //
582            // NOTE: we assume that P[i,j] equals to 0 or 1,
583            //       intermediate values are not allowed.
584            //
585            if( tiecount==1 )
586            {
587                info = -3;
588                return;
589            }
590           
591            //
592            // General case, number of ties > 1
593            //
594            // NOTE: we assume that P[i,j] equals to 0 or 1,
595            //       intermediate values are not allowed.
596            //
597            pal = 0;
598            pbl = 0;
599            par = 0;
600            pbr = 0;
601            for(i=0; i<=n-1; i++)
602            {
603                if( c[i]==0 )
604                {
605                    par = par+1;
606                }
607                if( c[i]==1 )
608                {
609                    pbr = pbr+1;
610                }
611            }
612            koptimal = -1;
613            cvoptimal = AP.Math.MaxRealNumber;
614            for(k=0; k<=tiecount-2; k++)
615            {
616               
617                //
618                // first, obtain information about K-th tie which is
619                // moved from R-part to L-part
620                //
621                pak = 0;
622                pbk = 0;
623                for(i=ties[k]; i<=ties[k+1]-1; i++)
624                {
625                    if( c[i]==0 )
626                    {
627                        pak = pak+1;
628                    }
629                    if( c[i]==1 )
630                    {
631                        pbk = pbk+1;
632                    }
633                }
634               
635                //
636                // Calculate cross-validation CE
637                //
638                cv = 0;
639                cv = cv-xlny(pal+pak, (pal+pak)/(pal+pak+pbl+pbk+1));
640                cv = cv-xlny(pbl+pbk, (pbl+pbk)/(pal+pak+1+pbl+pbk));
641                cv = cv-xlny(par-pak, (par-pak)/(par-pak+pbr-pbk+1));
642                cv = cv-xlny(pbr-pbk, (pbr-pbk)/(par-pak+1+pbr-pbk));
643               
644                //
645                // Compare with best
646                //
647                if( (double)(cv)<(double)(cvoptimal) )
648                {
649                    cvoptimal = cv;
650                    koptimal = k;
651                }
652               
653                //
654                // update
655                //
656                pal = pal+pak;
657                pbl = pbl+pbk;
658                par = par-pak;
659                pbr = pbr-pbk;
660            }
661            cve = cvoptimal;
662            threshold = 0.5*(a[ties[koptimal]]+a[ties[koptimal+1]]);
663            pal = 0;
664            pbl = 0;
665            par = 0;
666            pbr = 0;
667            for(i=0; i<=n-1; i++)
668            {
669                if( (double)(a[i])<(double)(threshold) )
670                {
671                    if( c[i]==0 )
672                    {
673                        pal = pal+1;
674                    }
675                    else
676                    {
677                        pbl = pbl+1;
678                    }
679                }
680                else
681                {
682                    if( c[i]==0 )
683                    {
684                        par = par+1;
685                    }
686                    else
687                    {
688                        pbr = pbr+1;
689                    }
690                }
691            }
692            s = pal+pbl;
693            pal = pal/s;
694            pbl = pbl/s;
695            s = par+pbr;
696            par = par/s;
697            pbr = pbr/s;
698        }
699
700
701        /*************************************************************************
702        Optimal partition, internal subroutine. Fast version.
703
704        Accepts:
705            A       array[0..N-1]       array of attributes     array[0..N-1]
706            C       array[0..N-1]       array of class labels
707            TiesBuf array[0..N]         temporaries (ties)
708            CntBuf  array[0..2*NC-1]    temporaries (counts)
709            Alpha                       centering factor (0<=alpha<=1, recommended value - 0.05)
710           
711        Output:
712            Info    error code (">0"=OK, "<0"=bad)
713            RMS     training set RMS error
714            CVRMS   leave-one-out RMS error
715           
716        Note:
717            content of all arrays is changed by subroutine
718
719          -- ALGLIB --
720             Copyright 11.12.2008 by Bochkanov Sergey
721        *************************************************************************/
722        public static void dsoptimalsplit2fast(ref double[] a,
723            ref int[] c,
724            ref int[] tiesbuf,
725            ref int[] cntbuf,
726            int n,
727            int nc,
728            double alpha,
729            ref int info,
730            ref double threshold,
731            ref double rms,
732            ref double cvrms)
733        {
734            int i = 0;
735            int k = 0;
736            int cl = 0;
737            int tiecount = 0;
738            double cbest = 0;
739            double cc = 0;
740            int koptimal = 0;
741            int sl = 0;
742            int sr = 0;
743            double v = 0;
744            double w = 0;
745            double x = 0;
746
747           
748            //
749            // Test for errors in inputs
750            //
751            if( n<=0 | nc<2 )
752            {
753                info = -1;
754                return;
755            }
756            for(i=0; i<=n-1; i++)
757            {
758                if( c[i]<0 | c[i]>=nc )
759                {
760                    info = -2;
761                    return;
762                }
763            }
764            info = 1;
765           
766            //
767            // Tie
768            //
769            dstiefasti(ref a, ref c, n, ref tiesbuf, ref tiecount);
770           
771            //
772            // Special case: number of ties is 1.
773            //
774            if( tiecount==1 )
775            {
776                info = -3;
777                return;
778            }
779           
780            //
781            // General case, number of ties > 1
782            //
783            for(i=0; i<=2*nc-1; i++)
784            {
785                cntbuf[i] = 0;
786            }
787            for(i=0; i<=n-1; i++)
788            {
789                cntbuf[nc+c[i]] = cntbuf[nc+c[i]]+1;
790            }
791            koptimal = -1;
792            threshold = a[n-1];
793            cbest = AP.Math.MaxRealNumber;
794            sl = 0;
795            sr = n;
796            for(k=0; k<=tiecount-2; k++)
797            {
798               
799                //
800                // first, move Kth tie from right to left
801                //
802                for(i=tiesbuf[k]; i<=tiesbuf[k+1]-1; i++)
803                {
804                    cl = c[i];
805                    cntbuf[cl] = cntbuf[cl]+1;
806                    cntbuf[nc+cl] = cntbuf[nc+cl]-1;
807                }
808                sl = sl+(tiesbuf[k+1]-tiesbuf[k]);
809                sr = sr-(tiesbuf[k+1]-tiesbuf[k]);
810               
811                //
812                // Calculate RMS error
813                //
814                v = 0;
815                for(i=0; i<=nc-1; i++)
816                {
817                    w = cntbuf[i];
818                    v = v+w*AP.Math.Sqr(w/sl-1);
819                    v = v+(sl-w)*AP.Math.Sqr(w/sl);
820                    w = cntbuf[nc+i];
821                    v = v+w*AP.Math.Sqr(w/sr-1);
822                    v = v+(sr-w)*AP.Math.Sqr(w/sr);
823                }
824                v = Math.Sqrt(v/(nc*n));
825               
826                //
827                // Compare with best
828                //
829                x = (double)(2*sl)/((double)(sl+sr))-1;
830                cc = v*(1-alpha+alpha*AP.Math.Sqr(x));
831                if( (double)(cc)<(double)(cbest) )
832                {
833                   
834                    //
835                    // store split
836                    //
837                    rms = v;
838                    koptimal = k;
839                    cbest = cc;
840                   
841                    //
842                    // calculate CVRMS error
843                    //
844                    cvrms = 0;
845                    for(i=0; i<=nc-1; i++)
846                    {
847                        if( sl>1 )
848                        {
849                            w = cntbuf[i];
850                            cvrms = cvrms+w*AP.Math.Sqr((w-1)/(sl-1)-1);
851                            cvrms = cvrms+(sl-w)*AP.Math.Sqr(w/(sl-1));
852                        }
853                        else
854                        {
855                            w = cntbuf[i];
856                            cvrms = cvrms+w*AP.Math.Sqr((double)(1)/(double)(nc)-1);
857                            cvrms = cvrms+(sl-w)*AP.Math.Sqr((double)(1)/(double)(nc));
858                        }
859                        if( sr>1 )
860                        {
861                            w = cntbuf[nc+i];
862                            cvrms = cvrms+w*AP.Math.Sqr((w-1)/(sr-1)-1);
863                            cvrms = cvrms+(sr-w)*AP.Math.Sqr(w/(sr-1));
864                        }
865                        else
866                        {
867                            w = cntbuf[nc+i];
868                            cvrms = cvrms+w*AP.Math.Sqr((double)(1)/(double)(nc)-1);
869                            cvrms = cvrms+(sr-w)*AP.Math.Sqr((double)(1)/(double)(nc));
870                        }
871                    }
872                    cvrms = Math.Sqrt(cvrms/(nc*n));
873                }
874            }
875           
876            //
877            // Calculate threshold.
878            // Code is a bit complicated because there can be such
879            // numbers that 0.5(A+B) equals to A or B (if A-B=epsilon)
880            //
881            threshold = 0.5*(a[tiesbuf[koptimal]]+a[tiesbuf[koptimal+1]]);
882            if( (double)(threshold)<=(double)(a[tiesbuf[koptimal]]) )
883            {
884                threshold = a[tiesbuf[koptimal+1]];
885            }
886        }
887
888
889        /*************************************************************************
890        Automatic non-optimal discretization, internal subroutine.
891
892          -- ALGLIB --
893             Copyright 22.05.2008 by Bochkanov Sergey
894        *************************************************************************/
895        public static void dssplitk(double[] a,
896            int[] c,
897            int n,
898            int nc,
899            int kmax,
900            ref int info,
901            ref double[] thresholds,
902            ref int ni,
903            ref double cve)
904        {
905            int i = 0;
906            int j = 0;
907            int j1 = 0;
908            int k = 0;
909            int[] ties = new int[0];
910            int tiecount = 0;
911            int[] p1 = new int[0];
912            int[] p2 = new int[0];
913            int[] cnt = new int[0];
914            double v2 = 0;
915            int bestk = 0;
916            double bestcve = 0;
917            int[] bestsizes = new int[0];
918            double curcve = 0;
919            int[] cursizes = new int[0];
920
921            a = (double[])a.Clone();
922            c = (int[])c.Clone();
923
924           
925            //
926            // Test for errors in inputs
927            //
928            if( n<=0 | nc<2 | kmax<2 )
929            {
930                info = -1;
931                return;
932            }
933            for(i=0; i<=n-1; i++)
934            {
935                if( c[i]<0 | c[i]>=nc )
936                {
937                    info = -2;
938                    return;
939                }
940            }
941            info = 1;
942           
943            //
944            // Tie
945            //
946            dstie(ref a, n, ref ties, ref tiecount, ref p1, ref p2);
947            for(i=0; i<=n-1; i++)
948            {
949                if( p2[i]!=i )
950                {
951                    k = c[i];
952                    c[i] = c[p2[i]];
953                    c[p2[i]] = k;
954                }
955            }
956           
957            //
958            // Special cases
959            //
960            if( tiecount==1 )
961            {
962                info = -3;
963                return;
964            }
965           
966            //
967            // General case:
968            // 0. allocate arrays
969            //
970            kmax = Math.Min(kmax, tiecount);
971            bestsizes = new int[kmax-1+1];
972            cursizes = new int[kmax-1+1];
973            cnt = new int[nc-1+1];
974           
975            //
976            // General case:
977            // 1. prepare "weak" solution (two subintervals, divided at median)
978            //
979            v2 = AP.Math.MaxRealNumber;
980            j = -1;
981            for(i=1; i<=tiecount-1; i++)
982            {
983                if( (double)(Math.Abs(ties[i]-0.5*(n-1)))<(double)(v2) )
984                {
985                    v2 = Math.Abs(ties[i]-0.5*n);
986                    j = i;
987                }
988            }
989            System.Diagnostics.Debug.Assert(j>0, "DSSplitK: internal error #1!");
990            bestk = 2;
991            bestsizes[0] = ties[j];
992            bestsizes[1] = n-j;
993            bestcve = 0;
994            for(i=0; i<=nc-1; i++)
995            {
996                cnt[i] = 0;
997            }
998            for(i=0; i<=j-1; i++)
999            {
1000                tieaddc(ref c, ref ties, i, nc, ref cnt);
1001            }
1002            bestcve = bestcve+getcv(ref cnt, nc);
1003            for(i=0; i<=nc-1; i++)
1004            {
1005                cnt[i] = 0;
1006            }
1007            for(i=j; i<=tiecount-1; i++)
1008            {
1009                tieaddc(ref c, ref ties, i, nc, ref cnt);
1010            }
1011            bestcve = bestcve+getcv(ref cnt, nc);
1012           
1013            //
1014            // General case:
1015            // 2. Use greedy algorithm to find sub-optimal split in O(KMax*N) time
1016            //
1017            for(k=2; k<=kmax; k++)
1018            {
1019               
1020                //
1021                // Prepare greedy K-interval split
1022                //
1023                for(i=0; i<=k-1; i++)
1024                {
1025                    cursizes[i] = 0;
1026                }
1027                i = 0;
1028                j = 0;
1029                while( j<=tiecount-1 & i<=k-1 )
1030                {
1031                   
1032                    //
1033                    // Rule: I-th bin is empty, fill it
1034                    //
1035                    if( cursizes[i]==0 )
1036                    {
1037                        cursizes[i] = ties[j+1]-ties[j];
1038                        j = j+1;
1039                        continue;
1040                    }
1041                   
1042                    //
1043                    // Rule: (K-1-I) bins left, (K-1-I) ties left (1 tie per bin); next bin
1044                    //
1045                    if( tiecount-j==k-1-i )
1046                    {
1047                        i = i+1;
1048                        continue;
1049                    }
1050                   
1051                    //
1052                    // Rule: last bin, always place in current
1053                    //
1054                    if( i==k-1 )
1055                    {
1056                        cursizes[i] = cursizes[i]+ties[j+1]-ties[j];
1057                        j = j+1;
1058                        continue;
1059                    }
1060                   
1061                    //
1062                    // Place J-th tie in I-th bin, or leave for I+1-th bin.
1063                    //
1064                    if( (double)(Math.Abs(cursizes[i]+ties[j+1]-ties[j]-(double)(n)/(double)(k)))<(double)(Math.Abs(cursizes[i]-(double)(n)/(double)(k))) )
1065                    {
1066                        cursizes[i] = cursizes[i]+ties[j+1]-ties[j];
1067                        j = j+1;
1068                    }
1069                    else
1070                    {
1071                        i = i+1;
1072                    }
1073                }
1074                System.Diagnostics.Debug.Assert(cursizes[k-1]!=0 & j==tiecount, "DSSplitK: internal error #1");
1075               
1076                //
1077                // Calculate CVE
1078                //
1079                curcve = 0;
1080                j = 0;
1081                for(i=0; i<=k-1; i++)
1082                {
1083                    for(j1=0; j1<=nc-1; j1++)
1084                    {
1085                        cnt[j1] = 0;
1086                    }
1087                    for(j1=j; j1<=j+cursizes[i]-1; j1++)
1088                    {
1089                        cnt[c[j1]] = cnt[c[j1]]+1;
1090                    }
1091                    curcve = curcve+getcv(ref cnt, nc);
1092                    j = j+cursizes[i];
1093                }
1094               
1095                //
1096                // Choose best variant
1097                //
1098                if( (double)(curcve)<(double)(bestcve) )
1099                {
1100                    for(i=0; i<=k-1; i++)
1101                    {
1102                        bestsizes[i] = cursizes[i];
1103                    }
1104                    bestcve = curcve;
1105                    bestk = k;
1106                }
1107            }
1108           
1109            //
1110            // Transform from sizes to thresholds
1111            //
1112            cve = bestcve;
1113            ni = bestk;
1114            thresholds = new double[ni-2+1];
1115            j = bestsizes[0];
1116            for(i=1; i<=bestk-1; i++)
1117            {
1118                thresholds[i-1] = 0.5*(a[j-1]+a[j]);
1119                j = j+bestsizes[i];
1120            }
1121        }
1122
1123
1124        /*************************************************************************
1125        Automatic optimal discretization, internal subroutine.
1126
1127          -- ALGLIB --
1128             Copyright 22.05.2008 by Bochkanov Sergey
1129        *************************************************************************/
1130        public static void dsoptimalsplitk(double[] a,
1131            int[] c,
1132            int n,
1133            int nc,
1134            int kmax,
1135            ref int info,
1136            ref double[] thresholds,
1137            ref int ni,
1138            ref double cve)
1139        {
1140            int i = 0;
1141            int j = 0;
1142            int s = 0;
1143            int jl = 0;
1144            int jr = 0;
1145            double v1 = 0;
1146            double v2 = 0;
1147            double v3 = 0;
1148            double v4 = 0;
1149            int[] ties = new int[0];
1150            int tiecount = 0;
1151            int[] p1 = new int[0];
1152            int[] p2 = new int[0];
1153            double cvtemp = 0;
1154            int[] cnt = new int[0];
1155            int[] cnt2 = new int[0];
1156            double[,] cv = new double[0,0];
1157            int[,] splits = new int[0,0];
1158            int k = 0;
1159            int koptimal = 0;
1160            double cvoptimal = 0;
1161
1162            a = (double[])a.Clone();
1163            c = (int[])c.Clone();
1164
1165           
1166            //
1167            // Test for errors in inputs
1168            //
1169            if( n<=0 | nc<2 | kmax<2 )
1170            {
1171                info = -1;
1172                return;
1173            }
1174            for(i=0; i<=n-1; i++)
1175            {
1176                if( c[i]<0 | c[i]>=nc )
1177                {
1178                    info = -2;
1179                    return;
1180                }
1181            }
1182            info = 1;
1183           
1184            //
1185            // Tie
1186            //
1187            dstie(ref a, n, ref ties, ref tiecount, ref p1, ref p2);
1188            for(i=0; i<=n-1; i++)
1189            {
1190                if( p2[i]!=i )
1191                {
1192                    k = c[i];
1193                    c[i] = c[p2[i]];
1194                    c[p2[i]] = k;
1195                }
1196            }
1197           
1198            //
1199            // Special cases
1200            //
1201            if( tiecount==1 )
1202            {
1203                info = -3;
1204                return;
1205            }
1206           
1207            //
1208            // General case
1209            // Use dynamic programming to find best split in O(KMax*NC*TieCount^2) time
1210            //
1211            kmax = Math.Min(kmax, tiecount);
1212            cv = new double[kmax-1+1, tiecount-1+1];
1213            splits = new int[kmax-1+1, tiecount-1+1];
1214            cnt = new int[nc-1+1];
1215            cnt2 = new int[nc-1+1];
1216            for(j=0; j<=nc-1; j++)
1217            {
1218                cnt[j] = 0;
1219            }
1220            for(j=0; j<=tiecount-1; j++)
1221            {
1222                tieaddc(ref c, ref ties, j, nc, ref cnt);
1223                splits[0,j] = 0;
1224                cv[0,j] = getcv(ref cnt, nc);
1225            }
1226            for(k=1; k<=kmax-1; k++)
1227            {
1228                for(j=0; j<=nc-1; j++)
1229                {
1230                    cnt[j] = 0;
1231                }
1232               
1233                //
1234                // Subtask size J in [K..TieCount-1]:
1235                // optimal K-splitting on ties from 0-th to J-th.
1236                //
1237                for(j=k; j<=tiecount-1; j++)
1238                {
1239                   
1240                    //
1241                    // Update Cnt - let it contain classes of ties from K-th to J-th
1242                    //
1243                    tieaddc(ref c, ref ties, j, nc, ref cnt);
1244                   
1245                    //
1246                    // Search for optimal split point S in [K..J]
1247                    //
1248                    for(i=0; i<=nc-1; i++)
1249                    {
1250                        cnt2[i] = cnt[i];
1251                    }
1252                    cv[k,j] = cv[k-1,j-1]+getcv(ref cnt2, nc);
1253                    splits[k,j] = j;
1254                    for(s=k+1; s<=j; s++)
1255                    {
1256                       
1257                        //
1258                        // Update Cnt2 - let it contain classes of ties from S-th to J-th
1259                        //
1260                        tiesubc(ref c, ref ties, s-1, nc, ref cnt2);
1261                       
1262                        //
1263                        // Calculate CVE
1264                        //
1265                        cvtemp = cv[k-1,s-1]+getcv(ref cnt2, nc);
1266                        if( (double)(cvtemp)<(double)(cv[k,j]) )
1267                        {
1268                            cv[k,j] = cvtemp;
1269                            splits[k,j] = s;
1270                        }
1271                    }
1272                }
1273            }
1274           
1275            //
1276            // Choose best partition, output result
1277            //
1278            koptimal = -1;
1279            cvoptimal = AP.Math.MaxRealNumber;
1280            for(k=0; k<=kmax-1; k++)
1281            {
1282                if( (double)(cv[k,tiecount-1])<(double)(cvoptimal) )
1283                {
1284                    cvoptimal = cv[k,tiecount-1];
1285                    koptimal = k;
1286                }
1287            }
1288            System.Diagnostics.Debug.Assert(koptimal>=0, "DSOptimalSplitK: internal error #1!");
1289            if( koptimal==0 )
1290            {
1291               
1292                //
1293                // Special case: best partition is one big interval.
1294                // Even 2-partition is not better.
1295                // This is possible when dealing with "weak" predictor variables.
1296                //
1297                // Make binary split as close to the median as possible.
1298                //
1299                v2 = AP.Math.MaxRealNumber;
1300                j = -1;
1301                for(i=1; i<=tiecount-1; i++)
1302                {
1303                    if( (double)(Math.Abs(ties[i]-0.5*(n-1)))<(double)(v2) )
1304                    {
1305                        v2 = Math.Abs(ties[i]-0.5*(n-1));
1306                        j = i;
1307                    }
1308                }
1309                System.Diagnostics.Debug.Assert(j>0, "DSOptimalSplitK: internal error #2!");
1310                thresholds = new double[0+1];
1311                thresholds[0] = 0.5*(a[ties[j-1]]+a[ties[j]]);
1312                ni = 2;
1313                cve = 0;
1314                for(i=0; i<=nc-1; i++)
1315                {
1316                    cnt[i] = 0;
1317                }
1318                for(i=0; i<=j-1; i++)
1319                {
1320                    tieaddc(ref c, ref ties, i, nc, ref cnt);
1321                }
1322                cve = cve+getcv(ref cnt, nc);
1323                for(i=0; i<=nc-1; i++)
1324                {
1325                    cnt[i] = 0;
1326                }
1327                for(i=j; i<=tiecount-1; i++)
1328                {
1329                    tieaddc(ref c, ref ties, i, nc, ref cnt);
1330                }
1331                cve = cve+getcv(ref cnt, nc);
1332            }
1333            else
1334            {
1335               
1336                //
1337                // General case: 2 or more intervals
1338                //
1339                thresholds = new double[koptimal-1+1];
1340                ni = koptimal+1;
1341                cve = cv[koptimal,tiecount-1];
1342                jl = splits[koptimal,tiecount-1];
1343                jr = tiecount-1;
1344                for(k=koptimal; k>=1; k--)
1345                {
1346                    thresholds[k-1] = 0.5*(a[ties[jl-1]]+a[ties[jl]]);
1347                    jr = jl-1;
1348                    jl = splits[k-1,jl-1];
1349                }
1350            }
1351        }
1352
1353
1354        /*************************************************************************
1355        Subroutine prepares K-fold split of the training set.
1356
1357        NOTES:
1358            "NClasses>0" means that we have classification task.
1359            "NClasses<0" means regression task with -NClasses real outputs.
1360
1361          -- ALGLIB --
1362             Copyright 11.01.2009 by Bochkanov Sergey
1363        *************************************************************************/
1364        private static void dskfoldsplit(ref double[,] xy,
1365            int npoints,
1366            int nclasses,
1367            int foldscount,
1368            bool stratifiedsplits,
1369            ref int[] folds)
1370        {
1371            int i = 0;
1372            int j = 0;
1373            int k = 0;
1374
1375           
1376            //
1377            // test parameters
1378            //
1379            System.Diagnostics.Debug.Assert(npoints>0, "DSKFoldSplit: wrong NPoints!");
1380            System.Diagnostics.Debug.Assert(nclasses>1 | nclasses<0, "DSKFoldSplit: wrong NClasses!");
1381            System.Diagnostics.Debug.Assert(foldscount>=2 & foldscount<=npoints, "DSKFoldSplit: wrong FoldsCount!");
1382            System.Diagnostics.Debug.Assert(!stratifiedsplits, "DSKFoldSplit: stratified splits are not supported!");
1383           
1384            //
1385            // Folds
1386            //
1387            folds = new int[npoints-1+1];
1388            for(i=0; i<=npoints-1; i++)
1389            {
1390                folds[i] = i*foldscount/npoints;
1391            }
1392            for(i=0; i<=npoints-2; i++)
1393            {
1394                j = i+AP.Math.RandomInteger(npoints-i);
1395                if( j!=i )
1396                {
1397                    k = folds[i];
1398                    folds[i] = folds[j];
1399                    folds[j] = k;
1400                }
1401            }
1402        }
1403
1404
1405        /*************************************************************************
1406        Internal function
1407        *************************************************************************/
1408        private static double xlny(double x,
1409            double y)
1410        {
1411            double result = 0;
1412
1413            if( (double)(x)==(double)(0) )
1414            {
1415                result = 0;
1416            }
1417            else
1418            {
1419                result = x*Math.Log(y);
1420            }
1421            return result;
1422        }
1423
1424
1425        /*************************************************************************
1426        Internal function,
1427        returns number of samples of class I in Cnt[I]
1428        *************************************************************************/
1429        private static double getcv(ref int[] cnt,
1430            int nc)
1431        {
1432            double result = 0;
1433            int i = 0;
1434            double s = 0;
1435
1436            s = 0;
1437            for(i=0; i<=nc-1; i++)
1438            {
1439                s = s+cnt[i];
1440            }
1441            result = 0;
1442            for(i=0; i<=nc-1; i++)
1443            {
1444                result = result-xlny(cnt[i], cnt[i]/(s+nc-1));
1445            }
1446            return result;
1447        }
1448
1449
1450        /*************************************************************************
1451        Internal function, adds number of samples of class I in tie NTie to Cnt[I]
1452        *************************************************************************/
1453        private static void tieaddc(ref int[] c,
1454            ref int[] ties,
1455            int ntie,
1456            int nc,
1457            ref int[] cnt)
1458        {
1459            int i = 0;
1460
1461            for(i=ties[ntie]; i<=ties[ntie+1]-1; i++)
1462            {
1463                cnt[c[i]] = cnt[c[i]]+1;
1464            }
1465        }
1466
1467
1468        /*************************************************************************
1469        Internal function, subtracts number of samples of class I in tie NTie to Cnt[I]
1470        *************************************************************************/
1471        private static void tiesubc(ref int[] c,
1472            ref int[] ties,
1473            int ntie,
1474            int nc,
1475            ref int[] cnt)
1476        {
1477            int i = 0;
1478
1479            for(i=ties[ntie]; i<=ties[ntie+1]-1; i++)
1480            {
1481                cnt[c[i]] = cnt[c[i]]-1;
1482            }
1483        }
1484
1485
1486        /*************************************************************************
1487        Internal function,
1488        returns number of samples of class I in Cnt[I]
1489        *************************************************************************/
1490        private static void tiegetc(ref int[] c,
1491            ref int[] ties,
1492            int ntie,
1493            int nc,
1494            ref int[] cnt)
1495        {
1496            int i = 0;
1497
1498            for(i=0; i<=nc-1; i++)
1499            {
1500                cnt[i] = 0;
1501            }
1502            for(i=ties[ntie]; i<=ties[ntie+1]-1; i++)
1503            {
1504                cnt[c[i]] = cnt[c[i]]+1;
1505            }
1506        }
1507    }
1508}
Note: See TracBrowser for help on using the repository browser.