Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.LinearRegression/3.2/alglib/reflections.cs @ 2211

Last change on this file since 2211 was 2154, checked in by gkronber, 15 years ago

Added linear regression plugin. #697

File size: 16.7 KB
Line 
1/*************************************************************************
2Copyright (c) 1992-2007 The University of Tennessee.  All rights reserved.
3
4Contributors:
5    * Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to
6      pseudocode.
7
8See subroutines comments for additional copyrights.
9
10Redistribution and use in source and binary forms, with or without
11modification, are permitted provided that the following conditions are
12met:
13
14- Redistributions of source code must retain the above copyright
15  notice, this list of conditions and the following disclaimer.
16
17- Redistributions in binary form must reproduce the above copyright
18  notice, this list of conditions and the following disclaimer listed
19  in this license in the documentation and/or other materials
20  provided with the distribution.
21
22- Neither the name of the copyright holders nor the names of its
23  contributors may be used to endorse or promote products derived from
24  this software without specific prior written permission.
25
26THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
27"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
28LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
29A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
30OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
31SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
32LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
33DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
34THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
35(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
36OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37*************************************************************************/
38
39using System;
40
41class reflections
42{
43    /*************************************************************************
44    Generation of an elementary reflection transformation
45
46    The subroutine generates elementary reflection H of order N, so that, for
47    a given X, the following equality holds true:
48
49        ( X(1) )   ( Beta )
50    H * (  ..  ) = (  0   )
51        ( X(n) )   (  0   )
52
53    where
54                  ( V(1) )
55    H = 1 - Tau * (  ..  ) * ( V(1), ..., V(n) )
56                  ( V(n) )
57
58    where the first component of vector V equals 1.
59
60    Input parameters:
61        X   -   vector. Array whose index ranges within [1..N].
62        N   -   reflection order.
63
64    Output parameters:
65        X   -   components from 2 to N are replaced with vector V.
66                The first component is replaced with parameter Beta.
67        Tau -   scalar value Tau. If X is a null vector, Tau equals 0,
68                otherwise 1 <= Tau <= 2.
69
70    This subroutine is the modification of the DLARFG subroutines from
71    the LAPACK library. It has a similar functionality except for the
72    fact that it doesn’t handle errors when the intermediate results
73    cause an overflow.
74
75
76    MODIFICATIONS:
77        24.12.2005 sign(Alpha) was replaced with an analogous to the Fortran SIGN code.
78
79      -- LAPACK auxiliary routine (version 3.0) --
80         Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
81         Courant Institute, Argonne National Lab, and Rice University
82         September 30, 1994
83    *************************************************************************/
84    public static void generatereflection(ref double[] x,
85        int n,
86        ref double tau)
87    {
88        int j = 0;
89        double alpha = 0;
90        double xnorm = 0;
91        double v = 0;
92        double beta = 0;
93        double mx = 0;
94        int i_ = 0;
95
96       
97        //
98        // Executable Statements ..
99        //
100        if( n<=1 )
101        {
102            tau = 0;
103            return;
104        }
105       
106        //
107        // XNORM = DNRM2( N-1, X, INCX )
108        //
109        alpha = x[1];
110        mx = 0;
111        for(j=2; j<=n; j++)
112        {
113            mx = Math.Max(Math.Abs(x[j]), mx);
114        }
115        xnorm = 0;
116        if( mx!=0 )
117        {
118            for(j=2; j<=n; j++)
119            {
120                xnorm = xnorm+AP.Math.Sqr(x[j]/mx);
121            }
122            xnorm = Math.Sqrt(xnorm)*mx;
123        }
124        if( xnorm==0 )
125        {
126           
127            //
128            // H  =  I
129            //
130            tau = 0;
131            return;
132        }
133       
134        //
135        // general case
136        //
137        mx = Math.Max(Math.Abs(alpha), Math.Abs(xnorm));
138        beta = -(mx*Math.Sqrt(AP.Math.Sqr(alpha/mx)+AP.Math.Sqr(xnorm/mx)));
139        if( alpha<0 )
140        {
141            beta = -beta;
142        }
143        tau = (beta-alpha)/beta;
144        v = 1/(alpha-beta);
145        for(i_=2; i_<=n;i_++)
146        {
147            x[i_] = v*x[i_];
148        }
149        x[1] = beta;
150    }
151
152
153    /*************************************************************************
154    Application of an elementary reflection to a rectangular matrix of size MxN
155
156    The algorithm pre-multiplies the matrix by an elementary reflection transformation
157    which is given by column V and scalar Tau (see the description of the
158    GenerateReflection procedure). Not the whole matrix but only a part of it
159    is transformed (rows from M1 to M2, columns from N1 to N2). Only the elements
160    of this submatrix are changed.
161
162    Input parameters:
163        C       -   matrix to be transformed.
164        Tau     -   scalar defining the transformation.
165        V       -   column defining the transformation.
166                    Array whose index ranges within [1..M2-M1+1].
167        M1, M2  -   range of rows to be transformed.
168        N1, N2  -   range of columns to be transformed.
169        WORK    -   working array whose indexes goes from N1 to N2.
170
171    Output parameters:
172        C       -   the result of multiplying the input matrix C by the
173                    transformation matrix which is given by Tau and V.
174                    If N1>N2 or M1>M2, C is not modified.
175
176      -- LAPACK auxiliary routine (version 3.0) --
177         Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
178         Courant Institute, Argonne National Lab, and Rice University
179         September 30, 1994
180    *************************************************************************/
181    public static void applyreflectionfromtheleft(ref double[,] c,
182        double tau,
183        ref double[] v,
184        int m1,
185        int m2,
186        int n1,
187        int n2,
188        ref double[] work)
189    {
190        double t = 0;
191        int i = 0;
192        int vm = 0;
193        int i_ = 0;
194
195        if( tau==0 | n1>n2 | m1>m2 )
196        {
197            return;
198        }
199       
200        //
201        // w := C' * v
202        //
203        vm = m2-m1+1;
204        for(i=n1; i<=n2; i++)
205        {
206            work[i] = 0;
207        }
208        for(i=m1; i<=m2; i++)
209        {
210            t = v[i+1-m1];
211            for(i_=n1; i_<=n2;i_++)
212            {
213                work[i_] = work[i_] + t*c[i,i_];
214            }
215        }
216       
217        //
218        // C := C - tau * v * w'
219        //
220        for(i=m1; i<=m2; i++)
221        {
222            t = v[i-m1+1]*tau;
223            for(i_=n1; i_<=n2;i_++)
224            {
225                c[i,i_] = c[i,i_] - t*work[i_];
226            }
227        }
228    }
229
230
231    /*************************************************************************
232    Application of an elementary reflection to a rectangular matrix of size MxN
233
234    The algorithm post-multiplies the matrix by an elementary reflection transformation
235    which is given by column V and scalar Tau (see the description of the
236    GenerateReflection procedure). Not the whole matrix but only a part of it
237    is transformed (rows from M1 to M2, columns from N1 to N2). Only the
238    elements of this submatrix are changed.
239
240    Input parameters:
241        C       -   matrix to be transformed.
242        Tau     -   scalar defining the transformation.
243        V       -   column defining the transformation.
244                    Array whose index ranges within [1..N2-N1+1].
245        M1, M2  -   range of rows to be transformed.
246        N1, N2  -   range of columns to be transformed.
247        WORK    -   working array whose indexes goes from M1 to M2.
248
249    Output parameters:
250        C       -   the result of multiplying the input matrix C by the
251                    transformation matrix which is given by Tau and V.
252                    If N1>N2 or M1>M2, C is not modified.
253
254      -- LAPACK auxiliary routine (version 3.0) --
255         Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
256         Courant Institute, Argonne National Lab, and Rice University
257         September 30, 1994
258    *************************************************************************/
259    public static void applyreflectionfromtheright(ref double[,] c,
260        double tau,
261        ref double[] v,
262        int m1,
263        int m2,
264        int n1,
265        int n2,
266        ref double[] work)
267    {
268        double t = 0;
269        int i = 0;
270        int vm = 0;
271        int i_ = 0;
272        int i1_ = 0;
273
274        if( tau==0 | n1>n2 | m1>m2 )
275        {
276            return;
277        }
278       
279        //
280        // w := C * v
281        //
282        vm = n2-n1+1;
283        for(i=m1; i<=m2; i++)
284        {
285            i1_ = (1)-(n1);
286            t = 0.0;
287            for(i_=n1; i_<=n2;i_++)
288            {
289                t += c[i,i_]*v[i_+i1_];
290            }
291            work[i] = t;
292        }
293       
294        //
295        // C := C - w * v'
296        //
297        for(i=m1; i<=m2; i++)
298        {
299            t = work[i]*tau;
300            i1_ = (1) - (n1);
301            for(i_=n1; i_<=n2;i_++)
302            {
303                c[i,i_] = c[i,i_] - t*v[i_+i1_];
304            }
305        }
306    }
307
308
309    private static void testreflections()
310    {
311        int i = 0;
312        int j = 0;
313        int n = 0;
314        int m = 0;
315        int maxmn = 0;
316        double[] x = new double[0];
317        double[] v = new double[0];
318        double[] work = new double[0];
319        double[,] h = new double[0,0];
320        double[,] a = new double[0,0];
321        double[,] b = new double[0,0];
322        double[,] c = new double[0,0];
323        double tmp = 0;
324        double beta = 0;
325        double tau = 0;
326        double err = 0;
327        double mer = 0;
328        double mel = 0;
329        double meg = 0;
330        int pass = 0;
331        int passcount = 0;
332        int i_ = 0;
333
334        passcount = 1000;
335        mer = 0;
336        mel = 0;
337        meg = 0;
338        for(pass=1; pass<=passcount; pass++)
339        {
340           
341            //
342            // Task
343            //
344            n = 1+AP.Math.RandomInteger(10);
345            m = 1+AP.Math.RandomInteger(10);
346            maxmn = Math.Max(m, n);
347           
348            //
349            // Initialize
350            //
351            x = new double[maxmn+1];
352            v = new double[maxmn+1];
353            work = new double[maxmn+1];
354            h = new double[maxmn+1, maxmn+1];
355            a = new double[maxmn+1, maxmn+1];
356            b = new double[maxmn+1, maxmn+1];
357            c = new double[maxmn+1, maxmn+1];
358           
359            //
360            // GenerateReflection
361            //
362            for(i=1; i<=n; i++)
363            {
364                x[i] = 2*AP.Math.RandomReal()-1;
365                v[i] = x[i];
366            }
367            generatereflection(ref v, n, ref tau);
368            beta = v[1];
369            v[1] = 1;
370            for(i=1; i<=n; i++)
371            {
372                for(j=1; j<=n; j++)
373                {
374                    if( i==j )
375                    {
376                        h[i,j] = 1-tau*v[i]*v[j];
377                    }
378                    else
379                    {
380                        h[i,j] = -(tau*v[i]*v[j]);
381                    }
382                }
383            }
384            err = 0;
385            for(i=1; i<=n; i++)
386            {
387                tmp = 0.0;
388                for(i_=1; i_<=n;i_++)
389                {
390                    tmp += h[i,i_]*x[i_];
391                }
392                if( i==1 )
393                {
394                    err = Math.Max(err, Math.Abs(tmp-beta));
395                }
396                else
397                {
398                    err = Math.Max(err, Math.Abs(tmp));
399                }
400            }
401            meg = Math.Max(meg, err);
402           
403            //
404            // ApplyReflectionFromTheLeft
405            //
406            for(i=1; i<=m; i++)
407            {
408                x[i] = 2*AP.Math.RandomReal()-1;
409                v[i] = x[i];
410            }
411            for(i=1; i<=m; i++)
412            {
413                for(j=1; j<=n; j++)
414                {
415                    a[i,j] = 2*AP.Math.RandomReal()-1;
416                    b[i,j] = a[i,j];
417                }
418            }
419            generatereflection(ref v, m, ref tau);
420            beta = v[1];
421            v[1] = 1;
422            applyreflectionfromtheleft(ref b, tau, ref v, 1, m, 1, n, ref work);
423            for(i=1; i<=m; i++)
424            {
425                for(j=1; j<=m; j++)
426                {
427                    if( i==j )
428                    {
429                        h[i,j] = 1-tau*v[i]*v[j];
430                    }
431                    else
432                    {
433                        h[i,j] = -(tau*v[i]*v[j]);
434                    }
435                }
436            }
437            for(i=1; i<=m; i++)
438            {
439                for(j=1; j<=n; j++)
440                {
441                    tmp = 0.0;
442                    for(i_=1; i_<=m;i_++)
443                    {
444                        tmp += h[i,i_]*a[i_,j];
445                    }
446                    c[i,j] = tmp;
447                }
448            }
449            err = 0;
450            for(i=1; i<=m; i++)
451            {
452                for(j=1; j<=n; j++)
453                {
454                    err = Math.Max(err, Math.Abs(b[i,j]-c[i,j]));
455                }
456            }
457            mel = Math.Max(mel, err);
458           
459            //
460            // ApplyReflectionFromTheRight
461            //
462            for(i=1; i<=n; i++)
463            {
464                x[i] = 2*AP.Math.RandomReal()-1;
465                v[i] = x[i];
466            }
467            for(i=1; i<=m; i++)
468            {
469                for(j=1; j<=n; j++)
470                {
471                    a[i,j] = 2*AP.Math.RandomReal()-1;
472                    b[i,j] = a[i,j];
473                }
474            }
475            generatereflection(ref v, n, ref tau);
476            beta = v[1];
477            v[1] = 1;
478            applyreflectionfromtheright(ref b, tau, ref v, 1, m, 1, n, ref work);
479            for(i=1; i<=n; i++)
480            {
481                for(j=1; j<=n; j++)
482                {
483                    if( i==j )
484                    {
485                        h[i,j] = 1-tau*v[i]*v[j];
486                    }
487                    else
488                    {
489                        h[i,j] = -(tau*v[i]*v[j]);
490                    }
491                }
492            }
493            for(i=1; i<=m; i++)
494            {
495                for(j=1; j<=n; j++)
496                {
497                    tmp = 0.0;
498                    for(i_=1; i_<=n;i_++)
499                    {
500                        tmp += a[i,i_]*h[i_,j];
501                    }
502                    c[i,j] = tmp;
503                }
504            }
505            err = 0;
506            for(i=1; i<=m; i++)
507            {
508                for(j=1; j<=n; j++)
509                {
510                    err = Math.Max(err, Math.Abs(b[i,j]-c[i,j]));
511                }
512            }
513            mer = Math.Max(mer, err);
514        }
515       
516        //
517        // Overflow crash test
518        //
519        x = new double[10+1];
520        v = new double[10+1];
521        for(i=1; i<=10; i++)
522        {
523            v[i] = AP.Math.MaxRealNumber*0.01*(2*AP.Math.RandomReal()-1);
524        }
525        generatereflection(ref v, 10, ref tau);
526        System.Console.Write("TESTING REFLECTIONS");
527        System.Console.WriteLine();
528        System.Console.Write("Pass count is ");
529        System.Console.Write("{0,0:d}",passcount);
530        System.Console.WriteLine();
531        System.Console.Write("Generate     absolute error is       ");
532        System.Console.Write("{0,5:E3}",meg);
533        System.Console.WriteLine();
534        System.Console.Write("Apply(Left)  absolute error is       ");
535        System.Console.Write("{0,5:E3}",mel);
536        System.Console.WriteLine();
537        System.Console.Write("Apply(Right) absolute error is       ");
538        System.Console.Write("{0,5:E3}",mer);
539        System.Console.WriteLine();
540        System.Console.Write("Overflow crash test passed");
541        System.Console.WriteLine();
542    }
543}
Note: See TracBrowser for help on using the repository browser.