Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.5.0/ALGLIB-3.5.0/ap.cs @ 7673

Last change on this file since 7673 was 7673, checked in by mkommend, 11 years ago

#1808: Added new alglib release 3.5.0 to the external libraries.

File size: 46.5 KB
Line 
1/*************************************************************************
2AP library
3Copyright (c) 2003-2009 Sergey Bochkanov (ALGLIB project).
4
5>>> SOURCE LICENSE >>>
6This program is free software; you can redistribute it and/or modify
7it under the terms of the GNU General Public License as published by
8the Free Software Foundation (www.fsf.org); either version 2 of the
9License, or (at your option) any later version.
10
11This program is distributed in the hope that it will be useful,
12but WITHOUT ANY WARRANTY; without even the implied warranty of
13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14GNU General Public License for more details.
15
16A copy of the GNU General Public License is available at
17http://www.fsf.org/licensing/licenses
18>>> END OF LICENSE >>>
19*************************************************************************/
20using System;
21public partial class alglib
22{
23    /********************************************************************
24    Callback definitions for optimizers/fitters/solvers.
25   
26    Callbacks for unparameterized (general) functions:
27    * ndimensional_func         calculates f(arg), stores result to func
28    * ndimensional_grad         calculates func = f(arg),
29                                grad[i] = df(arg)/d(arg[i])
30    * ndimensional_hess         calculates func = f(arg),
31                                grad[i] = df(arg)/d(arg[i]),
32                                hess[i,j] = d2f(arg)/(d(arg[i])*d(arg[j]))
33   
34    Callbacks for systems of functions:
35    * ndimensional_fvec         calculates vector function f(arg),
36                                stores result to fi
37    * ndimensional_jac          calculates f[i] = fi(arg)
38                                jac[i,j] = df[i](arg)/d(arg[j])
39                               
40    Callbacks for  parameterized  functions,  i.e.  for  functions  which
41    depend on two vectors: P and Q.  Gradient  and Hessian are calculated
42    with respect to P only.
43    * ndimensional_pfunc        calculates f(p,q),
44                                stores result to func
45    * ndimensional_pgrad        calculates func = f(p,q),
46                                grad[i] = df(p,q)/d(p[i])
47    * ndimensional_phess        calculates func = f(p,q),
48                                grad[i] = df(p,q)/d(p[i]),
49                                hess[i,j] = d2f(p,q)/(d(p[i])*d(p[j]))
50
51    Callbacks for progress reports:
52    * ndimensional_rep          reports current position of optimization algo   
53   
54    Callbacks for ODE solvers:
55    * ndimensional_ode_rp       calculates dy/dx for given y[] and x
56   
57    Callbacks for integrators:
58    * integrator1_func          calculates f(x) for given x
59                                (additional parameters xminusa and bminusx
60                                contain x-a and b-x)
61    ********************************************************************/
62    public delegate void ndimensional_func (double[] arg, ref double func, object obj);
63    public delegate void ndimensional_grad (double[] arg, ref double func, double[] grad, object obj);
64    public delegate void ndimensional_hess (double[] arg, ref double func, double[] grad, double[,] hess, object obj);
65   
66    public delegate void ndimensional_fvec (double[] arg, double[] fi, object obj);
67    public delegate void ndimensional_jac  (double[] arg, double[] fi, double[,] jac, object obj);
68   
69    public delegate void ndimensional_pfunc(double[] p, double[] q, ref double func, object obj);
70    public delegate void ndimensional_pgrad(double[] p, double[] q, ref double func, double[] grad, object obj);
71    public delegate void ndimensional_phess(double[] p, double[] q, ref double func, double[] grad, double[,] hess, object obj);
72   
73    public delegate void ndimensional_rep(double[] arg, double func, object obj);
74
75    public delegate void ndimensional_ode_rp (double[] y, double x, double[] dy, object obj);
76
77    public delegate void integrator1_func (double x, double xminusa, double bminusx, ref double f, object obj);
78
79    /********************************************************************
80    Class defining a complex number with double precision.
81    ********************************************************************/
82    public struct complex
83    {
84        public double x;
85        public double y;
86
87        public complex(double _x)
88        {
89            x = _x;
90            y = 0;
91        }
92        public complex(double _x, double _y)
93        {
94            x = _x;
95            y = _y;
96        }
97        public static implicit operator complex(double _x)
98        {
99            return new complex(_x);
100        }
101        public static bool operator==(complex lhs, complex rhs)
102        {
103            return ((double)lhs.x==(double)rhs.x) & ((double)lhs.y==(double)rhs.y);
104        }
105        public static bool operator!=(complex lhs, complex rhs)
106        {
107            return ((double)lhs.x!=(double)rhs.x) | ((double)lhs.y!=(double)rhs.y);
108        }
109        public static complex operator+(complex lhs)
110        {
111            return lhs;
112        }
113        public static complex operator-(complex lhs)
114        {
115            return new complex(-lhs.x,-lhs.y);
116        }
117        public static complex operator+(complex lhs, complex rhs)
118        {
119            return new complex(lhs.x+rhs.x,lhs.y+rhs.y);
120        }
121        public static complex operator-(complex lhs, complex rhs)
122        {
123            return new complex(lhs.x-rhs.x,lhs.y-rhs.y);
124        }
125        public static complex operator*(complex lhs, complex rhs)
126        {
127            return new complex(lhs.x*rhs.x-lhs.y*rhs.y, lhs.x*rhs.y+lhs.y*rhs.x);
128        }
129        public static complex operator/(complex lhs, complex rhs)
130        {
131            complex result;
132            double e;
133            double f;
134            if( System.Math.Abs(rhs.y)<System.Math.Abs(rhs.x) )
135            {
136                e = rhs.y/rhs.x;
137                f = rhs.x+rhs.y*e;
138                result.x = (lhs.x+lhs.y*e)/f;
139                result.y = (lhs.y-lhs.x*e)/f;
140            }
141            else
142            {
143                e = rhs.x/rhs.y;
144                f = rhs.y+rhs.x*e;
145                result.x = (lhs.y+lhs.x*e)/f;
146                result.y = (-lhs.x+lhs.y*e)/f;
147            }
148            return result;
149        }
150        public override int GetHashCode()
151        {
152            return x.GetHashCode() ^ y.GetHashCode();
153        }
154        public override bool Equals(object obj)
155        {
156            if( obj is byte)
157                return Equals(new complex((byte)obj));
158            if( obj is sbyte)
159                return Equals(new complex((sbyte)obj));
160            if( obj is short)
161                return Equals(new complex((short)obj));
162            if( obj is ushort)
163                return Equals(new complex((ushort)obj));
164            if( obj is int)
165                return Equals(new complex((int)obj));
166            if( obj is uint)
167                return Equals(new complex((uint)obj));
168            if( obj is long)
169                return Equals(new complex((long)obj));
170            if( obj is ulong)
171                return Equals(new complex((ulong)obj));
172            if( obj is float)
173                return Equals(new complex((float)obj));
174            if( obj is double)
175                return Equals(new complex((double)obj));
176            if( obj is decimal)
177                return Equals(new complex((double)(decimal)obj));
178            return base.Equals(obj);
179        }   
180    }   
181   
182    /********************************************************************
183    Class defining an ALGLIB exception
184    ********************************************************************/
185    public class alglibexception : System.Exception
186    {
187        public string msg;
188        public alglibexception(string s)
189        {
190            msg = s;
191        }
192       
193    }
194   
195    /********************************************************************
196    reverse communication structure
197    ********************************************************************/
198    public class rcommstate
199    {
200        public rcommstate()
201        {
202            stage = -1;
203            ia = new int[0];
204            ba = new bool[0];
205            ra = new double[0];
206            ca = new alglib.complex[0];
207        }
208        public int stage;
209        public int[] ia;
210        public bool[] ba;
211        public double[] ra;
212        public alglib.complex[] ca;
213    };
214
215    /********************************************************************
216    internal functions
217    ********************************************************************/
218    public class ap
219    {
220        public static int len<T>(T[] a)
221        { return a.Length; }
222        public static int rows<T>(T[,] a)
223        { return a.GetLength(0); }
224        public static int cols<T>(T[,] a)
225        { return a.GetLength(1); }
226        public static void swap<T>(ref T a, ref T b)
227        {
228            T t = a;
229            a = b;
230            b = t;
231        }
232       
233        public static void assert(bool cond, string s)
234        {
235            if( !cond )
236                throw new alglibexception(s);
237        }
238       
239        public static void assert(bool cond)
240        {
241            assert(cond, "ALGLIB: assertion failed");
242        }
243       
244        /****************************************************************
245        returns dps (digits-of-precision) value corresponding to threshold.
246        dps(0.9)  = dps(0.5)  = dps(0.1) = 0
247        dps(0.09) = dps(0.05) = dps(0.01) = 1
248        and so on
249        ****************************************************************/
250        public static int threshold2dps(double threshold)
251        {
252            int result = 0;
253            double t;
254            for (result = 0, t = 1; t / 10 > threshold*(1+1E-10); result++, t /= 10) ;
255            return result;
256        }
257
258        /****************************************************************
259        prints formatted complex
260        ****************************************************************/
261        public static string format(complex a, int _dps)
262        {
263            int dps = Math.Abs(_dps);
264            string fmt = _dps>=0 ? "F" : "E";
265            string fmtx = String.Format("{{0:"+fmt+"{0}}}", dps);
266            string fmty = String.Format("{{0:"+fmt+"{0}}}", dps);
267            string result = String.Format(fmtx, a.x) + (a.y >= 0 ? "+" : "-") + String.Format(fmty, Math.Abs(a.y)) + "i";
268            result = result.Replace(',', '.');
269            return result;
270        }
271
272        /****************************************************************
273        prints formatted array
274        ****************************************************************/
275        public static string format(bool[] a)
276        {
277            string[] result = new string[len(a)];
278            int i;
279            for(i=0; i<len(a); i++)
280                if( a[i] )
281                    result[i] = "true";
282                else
283                    result[i] = "false";
284            return "{"+String.Join(",",result)+"}";
285        }
286       
287        /****************************************************************
288        prints formatted array
289        ****************************************************************/
290        public static string format(int[] a)
291        {
292            string[] result = new string[len(a)];
293            int i;
294            for (i = 0; i < len(a); i++)
295                result[i] = a[i].ToString();
296            return "{" + String.Join(",", result) + "}";
297        }
298
299        /****************************************************************
300        prints formatted array
301        ****************************************************************/
302        public static string format(double[] a, int _dps)
303        {
304            int dps = Math.Abs(_dps);
305            string sfmt = _dps >= 0 ? "F" : "E";
306            string fmt = String.Format("{{0:" + sfmt + "{0}}}", dps);
307            string[] result = new string[len(a)];
308            int i;
309            for (i = 0; i < len(a); i++)
310            {
311                result[i] = String.Format(fmt, a[i]);
312                result[i] = result[i].Replace(',', '.');
313            }
314            return "{" + String.Join(",", result) + "}";
315        }
316
317        /****************************************************************
318        prints formatted array
319        ****************************************************************/
320        public static string format(complex[] a, int _dps)
321        {
322            int dps = Math.Abs(_dps);
323            string fmt = _dps >= 0 ? "F" : "E";
324            string fmtx = String.Format("{{0:"+fmt+"{0}}}", dps);
325            string fmty = String.Format("{{0:"+fmt+"{0}}}", dps);
326            string[] result = new string[len(a)];
327            int i;
328            for (i = 0; i < len(a); i++)
329            {
330                result[i] = String.Format(fmtx, a[i].x) + (a[i].y >= 0 ? "+" : "-") + String.Format(fmty, Math.Abs(a[i].y)) + "i";
331                result[i] = result[i].Replace(',', '.');
332            }
333            return "{" + String.Join(",", result) + "}";
334        }
335
336        /****************************************************************
337        prints formatted matrix
338        ****************************************************************/
339        public static string format(bool[,] a)
340        {
341            int i, j, m, n;
342            n = cols(a);
343            m = rows(a);
344            bool[] line = new bool[n];
345            string[] result = new string[m];
346            for (i = 0; i < m; i++)
347            {
348                for (j = 0; j < n; j++)
349                    line[j] = a[i, j];
350                result[i] = format(line);
351            }
352            return "{" + String.Join(",", result) + "}";
353        }
354
355        /****************************************************************
356        prints formatted matrix
357        ****************************************************************/
358        public static string format(int[,] a)
359        {
360            int i, j, m, n;
361            n = cols(a);
362            m = rows(a);
363            int[] line = new int[n];
364            string[] result = new string[m];
365            for (i = 0; i < m; i++)
366            {
367                for (j = 0; j < n; j++)
368                    line[j] = a[i, j];
369                result[i] = format(line);
370            }
371            return "{" + String.Join(",", result) + "}";
372        }
373
374        /****************************************************************
375        prints formatted matrix
376        ****************************************************************/
377        public static string format(double[,] a, int dps)
378        {
379            int i, j, m, n;
380            n = cols(a);
381            m = rows(a);
382            double[] line = new double[n];
383            string[] result = new string[m];
384            for (i = 0; i < m; i++)
385            {
386                for (j = 0; j < n; j++)
387                    line[j] = a[i, j];
388                result[i] = format(line, dps);
389            }
390            return "{" + String.Join(",", result) + "}";
391        }
392
393        /****************************************************************
394        prints formatted matrix
395        ****************************************************************/
396        public static string format(complex[,] a, int dps)
397        {
398            int i, j, m, n;
399            n = cols(a);
400            m = rows(a);
401            complex[] line = new complex[n];
402            string[] result = new string[m];
403            for (i = 0; i < m; i++)
404            {
405                for (j = 0; j < n; j++)
406                    line[j] = a[i, j];
407                result[i] = format(line, dps);
408            }
409            return "{" + String.Join(",", result) + "}";
410        }
411
412        /****************************************************************
413        checks that matrix is symmetric.
414        max|A-A^T| is calculated; if it is within 1.0E-14 of max|A|,
415        matrix is considered symmetric
416        ****************************************************************/
417        public static bool issymmetric(double[,] a)
418        {
419            int i, j, n;
420            double err, mx, v1, v2;
421            if( rows(a)!=cols(a) )
422                return false;
423            n = rows(a);
424            if( n==0 )
425                return true;
426            mx = 0;
427            err = 0;
428            for( i=0; i<n; i++)
429            {
430                for(j=i+1; j<n; j++)
431                {
432                    v1 = a[i,j];
433                    v2 = a[j,i];
434                    if( !math.isfinite(v1) )
435                        return false;
436                    if( !math.isfinite(v2) )
437                        return false;
438                    err = Math.Max(err, Math.Abs(v1-v2));
439                    mx  = Math.Max(mx,  Math.Abs(v1));
440                    mx  = Math.Max(mx,  Math.Abs(v2));
441                }
442                v1 = a[i,i];
443                if( !math.isfinite(v1) )
444                    return false;
445                mx = Math.Max(mx, Math.Abs(v1));
446            }
447            if( mx==0 )
448                return true;
449            return err/mx<=1.0E-14;
450        }
451       
452        /****************************************************************
453        checks that matrix is Hermitian.
454        max|A-A^H| is calculated; if it is within 1.0E-14 of max|A|,
455        matrix is considered Hermitian
456        ****************************************************************/
457        public static bool ishermitian(complex[,] a)
458        {
459            int i, j, n;
460            double err, mx;
461            complex v1, v2, vt;
462            if( rows(a)!=cols(a) )
463                return false;
464            n = rows(a);
465            if( n==0 )
466                return true;
467            mx = 0;
468            err = 0;
469            for( i=0; i<n; i++)
470            {
471                for(j=i+1; j<n; j++)
472                {
473                    v1 = a[i,j];
474                    v2 = a[j,i];
475                    if( !math.isfinite(v1.x) )
476                        return false;
477                    if( !math.isfinite(v1.y) )
478                        return false;
479                    if( !math.isfinite(v2.x) )
480                        return false;
481                    if( !math.isfinite(v2.y) )
482                        return false;
483                    vt.x = v1.x-v2.x;
484                    vt.y = v1.y+v2.y;
485                    err = Math.Max(err, math.abscomplex(vt));
486                    mx  = Math.Max(mx,  math.abscomplex(v1));
487                    mx  = Math.Max(mx,  math.abscomplex(v2));
488                }
489                v1 = a[i,i];
490                if( !math.isfinite(v1.x) )
491                    return false;
492                if( !math.isfinite(v1.y) )
493                    return false;
494                err = Math.Max(err, Math.Abs(v1.y));
495                mx = Math.Max(mx, math.abscomplex(v1));
496            }
497            if( mx==0 )
498                return true;
499            return err/mx<=1.0E-14;
500        }
501       
502       
503        /****************************************************************
504        Forces symmetricity by copying upper half of A to the lower one
505        ****************************************************************/
506        public static bool forcesymmetric(double[,] a)
507        {
508            int i, j, n;
509            if( rows(a)!=cols(a) )
510                return false;
511            n = rows(a);
512            if( n==0 )
513                return true;
514            for( i=0; i<n; i++)
515                for(j=i+1; j<n; j++)
516                    a[i,j] = a[j,i];
517            return true;
518        }
519       
520        /****************************************************************
521        Forces Hermiticity by copying upper half of A to the lower one
522        ****************************************************************/
523        public static bool forcehermitian(complex[,] a)
524        {
525            int i, j, n;
526            complex v;
527            if( rows(a)!=cols(a) )
528                return false;
529            n = rows(a);
530            if( n==0 )
531                return true;
532            for( i=0; i<n; i++)
533                for(j=i+1; j<n; j++)
534                {
535                    v = a[j,i];
536                    a[i,j].x = v.x;
537                    a[i,j].y = -v.y;
538                }
539            return true;
540        }
541    };
542   
543    /********************************************************************
544    math functions
545    ********************************************************************/
546    public class math
547    {
548        //public static System.Random RndObject = new System.Random(System.DateTime.Now.Millisecond);
549        public static System.Random rndobject = new System.Random(System.DateTime.Now.Millisecond + 1000*System.DateTime.Now.Second + 60*1000*System.DateTime.Now.Minute);
550
551        public const double machineepsilon = 5E-16;
552        public const double maxrealnumber = 1E300;
553        public const double minrealnumber = 1E-300;
554       
555        public static bool isfinite(double d)
556        {
557            return !System.Double.IsNaN(d) && !System.Double.IsInfinity(d);
558        }
559       
560        public static double randomreal()
561        {
562            double r = 0;
563            lock(rndobject){ r = rndobject.NextDouble(); }
564            return r;
565        }
566        public static int randominteger(int N)
567        {
568            int r = 0;
569            lock(rndobject){ r = rndobject.Next(N); }
570            return r;
571        }
572        public static double sqr(double X)
573        {
574            return X*X;
575        }       
576        public static double abscomplex(complex z)
577        {
578            double w;
579            double xabs;
580            double yabs;
581            double v;
582   
583            xabs = System.Math.Abs(z.x);
584            yabs = System.Math.Abs(z.y);
585            w = xabs>yabs ? xabs : yabs;
586            v = xabs<yabs ? xabs : yabs;
587            if( v==0 )
588                return w;
589            else
590            {
591                double t = v/w;
592                return w*System.Math.Sqrt(1+t*t);
593            }
594        }
595        public static complex conj(complex z)
596        {
597            return new complex(z.x, -z.y);
598        }   
599        public static complex csqr(complex z)
600        {
601            return new complex(z.x*z.x-z.y*z.y, 2*z.x*z.y);
602        }
603
604    }
605   
606    /********************************************************************
607    serializer object (should not be used directly)
608    ********************************************************************/
609    public class serializer
610    {
611        enum SMODE { DEFAULT, ALLOC, TO_STRING, FROM_STRING };
612        private const int SER_ENTRIES_PER_ROW = 5;
613        private const int SER_ENTRY_LENGTH    = 11;
614       
615        private SMODE mode;
616        private int entries_needed;
617        private int entries_saved;
618        private int bytes_asked;
619        private int bytes_written;
620        private int bytes_read;
621        private char[] out_str;
622        private char[] in_str;
623       
624        public serializer()
625        {
626            mode = SMODE.DEFAULT;
627            entries_needed = 0;
628            bytes_asked = 0;
629        }
630
631        public void alloc_start()
632        {
633            entries_needed = 0;
634            bytes_asked = 0;
635            mode = SMODE.ALLOC;
636        }
637
638        public void alloc_entry()
639        {
640            if( mode!=SMODE.ALLOC )
641                throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");
642            entries_needed++;
643        }
644
645        private int get_alloc_size()
646        {
647            int rows, lastrowsize, result;
648           
649            // check and change mode
650            if( mode!=SMODE.ALLOC )
651                throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");
652           
653            // if no entries needes (degenerate case)
654            if( entries_needed==0 )
655            {
656                bytes_asked = 1;
657                return bytes_asked;
658            }
659           
660            // non-degenerate case
661            rows = entries_needed/SER_ENTRIES_PER_ROW;
662            lastrowsize = SER_ENTRIES_PER_ROW;
663            if( entries_needed%SER_ENTRIES_PER_ROW!=0 )
664            {
665                lastrowsize = entries_needed%SER_ENTRIES_PER_ROW;
666                rows++;
667            }
668           
669            // calculate result size
670            result  = ((rows-1)*SER_ENTRIES_PER_ROW+lastrowsize)*SER_ENTRY_LENGTH;
671            result +=  (rows-1)*(SER_ENTRIES_PER_ROW-1)+(lastrowsize-1);
672            result += rows*2;
673            bytes_asked = result;
674            return result;
675        }
676
677        public void sstart_str()
678        {
679            int allocsize = get_alloc_size();
680           
681            // check and change mode
682            if( mode!=SMODE.ALLOC )
683                throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");
684            mode = SMODE.TO_STRING;
685           
686            // other preparations
687            out_str = new char[allocsize];
688            entries_saved = 0;
689            bytes_written = 0;
690        }
691
692        public void ustart_str(string s)
693        {
694            // check and change mode
695            if( mode!=SMODE.DEFAULT )
696                throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");
697            mode = SMODE.FROM_STRING;
698           
699            in_str = s.ToCharArray();
700            bytes_read = 0;
701        }
702
703        public void serialize_bool(bool v)
704        {
705            if( mode!=SMODE.TO_STRING )
706                throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");
707            bool2str(v, out_str, ref bytes_written);
708            entries_saved++;
709            if( entries_saved%SER_ENTRIES_PER_ROW!=0 )
710            {
711                out_str[bytes_written] = ' ';
712                bytes_written++;
713            }
714            else
715            {
716                out_str[bytes_written+0] = '\r';
717                out_str[bytes_written+1] = '\n';
718                bytes_written+=2;
719            }           
720        }
721
722        public void serialize_int(int v)
723        {
724            if( mode!=SMODE.TO_STRING )
725                throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");
726            int2str(v, out_str, ref bytes_written);
727            entries_saved++;
728            if( entries_saved%SER_ENTRIES_PER_ROW!=0 )
729            {
730                out_str[bytes_written] = ' ';
731                bytes_written++;
732            }
733            else
734            {
735                out_str[bytes_written+0] = '\r';
736                out_str[bytes_written+1] = '\n';
737                bytes_written+=2;
738            }
739        }
740
741        public void serialize_double(double v)
742        {
743            if( mode!=SMODE.TO_STRING )
744                throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");
745            double2str(v, out_str, ref bytes_written);
746            entries_saved++;
747            if( entries_saved%SER_ENTRIES_PER_ROW!=0 )
748            {
749                out_str[bytes_written] = ' ';
750                bytes_written++;
751            }
752            else
753            {
754                out_str[bytes_written+0] = '\r';
755                out_str[bytes_written+1] = '\n';
756                bytes_written+=2;
757            }
758        }
759
760        public bool unserialize_bool()
761        {
762            if( mode!=SMODE.FROM_STRING )
763                throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");
764            return str2bool(in_str, ref bytes_read);
765        }
766
767        public int unserialize_int()
768        {
769            if( mode!=SMODE.FROM_STRING )
770                throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");
771            return str2int(in_str, ref bytes_read);
772        }
773
774        public double unserialize_double()
775        {
776            if( mode!=SMODE.FROM_STRING )
777                throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");
778            return str2double(in_str, ref bytes_read);
779        }
780
781        public void stop()
782        {
783        }
784
785        public string get_string()
786        {
787            return new string(out_str, 0, bytes_written);
788        }
789
790
791        /************************************************************************
792        This function converts six-bit value (from 0 to 63)  to  character  (only
793        digits, lowercase and uppercase letters, minus and underscore are used).
794
795        If v is negative or greater than 63, this function returns '?'.
796        ************************************************************************/
797        private static char[] _sixbits2char_tbl = new char[64]{
798                '0', '1', '2', '3', '4', '5', '6', '7',
799                '8', '9', 'A', 'B', 'C', 'D', 'E', 'F',
800                'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N',
801                'O', 'P', 'Q', 'R', 'S', 'T', 'U', 'V',
802                'W', 'X', 'Y', 'Z', 'a', 'b', 'c', 'd',
803                'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l',
804                'm', 'n', 'o', 'p', 'q', 'r', 's', 't',
805                'u', 'v', 'w', 'x', 'y', 'z', '-', '_' };
806        private static char sixbits2char(int v)
807        {
808            if( v<0 || v>63 )
809                return '?';
810            return _sixbits2char_tbl[v];
811        }
812       
813        /************************************************************************
814        This function converts character to six-bit value (from 0 to 63).
815
816        This function is inverse of ae_sixbits2char()
817        If c is not correct character, this function returns -1.
818        ************************************************************************/
819        private static int[] _char2sixbits_tbl = new int[128] {
820            -1, -1, -1, -1, -1, -1, -1, -1,
821            -1, -1, -1, -1, -1, -1, -1, -1,
822            -1, -1, -1, -1, -1, -1, -1, -1,
823            -1, -1, -1, -1, -1, -1, -1, -1,
824            -1, -1, -1, -1, -1, -1, -1, -1,
825            -1, -1, -1, -1, -1, 62, -1, -1,
826             0,  1,  2,  3,  4,  5,  6,  7,
827             8,  9, -1, -1, -1, -1, -1, -1,
828            -1, 10, 11, 12, 13, 14, 15, 16,
829            17, 18, 19, 20, 21, 22, 23, 24,
830            25, 26, 27, 28, 29, 30, 31, 32,
831            33, 34, 35, -1, -1, -1, -1, 63,
832            -1, 36, 37, 38, 39, 40, 41, 42,
833            43, 44, 45, 46, 47, 48, 49, 50,
834            51, 52, 53, 54, 55, 56, 57, 58,
835            59, 60, 61, -1, -1, -1, -1, -1 };
836        private static int char2sixbits(char c)
837        {
838            return (c>=0 && c<127) ? _char2sixbits_tbl[c] : -1;
839        }
840       
841        /************************************************************************
842        This function converts three bytes (24 bits) to four six-bit values
843        (24 bits again).
844
845        src         array
846        src_offs    offset of three-bytes chunk
847        dst         array for ints
848        dst_offs    offset of four-ints chunk
849        ************************************************************************/
850        private static void threebytes2foursixbits(byte[] src, int src_offs, int[] dst, int dst_offs)
851        {
852            dst[dst_offs+0] =  src[src_offs+0] & 0x3F;
853            dst[dst_offs+1] = (src[src_offs+0]>>6) | ((src[src_offs+1]&0x0F)<<2);
854            dst[dst_offs+2] = (src[src_offs+1]>>4) | ((src[src_offs+2]&0x03)<<4);
855            dst[dst_offs+3] =  src[src_offs+2]>>2;
856        }
857
858        /************************************************************************
859        This function converts four six-bit values (24 bits) to three bytes
860        (24 bits again).
861
862        src         pointer to four ints
863        src_offs    offset of the chunk
864        dst         pointer to three bytes
865        dst_offs    offset of the chunk
866        ************************************************************************/
867        private static void foursixbits2threebytes(int[] src, int src_offs, byte[] dst, int dst_offs)
868        {
869            dst[dst_offs+0] =      (byte)(src[src_offs+0] | ((src[src_offs+1]&0x03)<<6));
870            dst[dst_offs+1] = (byte)((src[src_offs+1]>>2) | ((src[src_offs+2]&0x0F)<<4));
871            dst[dst_offs+2] = (byte)((src[src_offs+2]>>4) |  (src[src_offs+3]<<2));
872        }
873
874        /************************************************************************
875        This function serializes boolean value into buffer
876
877        v           boolean value to be serialized
878        buf         buffer, at least 11 characters wide
879        offs        offset in the buffer
880       
881        after return from this function, offs points to the char's past the value
882        being read.
883        ************************************************************************/
884        private static void bool2str(bool v, char[] buf, ref int offs)
885        {
886            char c = v ? '1' : '0';
887            int i;
888            for(i=0; i<SER_ENTRY_LENGTH; i++)
889                buf[offs+i] = c;
890            offs += SER_ENTRY_LENGTH;
891        }
892
893        /************************************************************************
894        This function unserializes boolean value from buffer
895
896        buf         buffer which contains value; leading spaces/tabs/newlines are
897                    ignored, traling spaces/tabs/newlines are treated as  end  of
898                    the boolean value.
899        offs        offset in the buffer
900       
901        after return from this function, offs points to the char's past the value
902        being read.
903
904        This function raises an error in case unexpected symbol is found
905        ************************************************************************/
906        private static bool str2bool(char[] buf, ref int offs)
907        {
908            bool was0, was1;
909            string emsg = "ALGLIB: unable to read boolean value from stream";
910           
911            was0 = false;
912            was1 = false;
913            while( buf[offs]==' ' || buf[offs]=='\t' || buf[offs]=='\n' || buf[offs]=='\r' )
914                offs++;
915            while( buf[offs]!=' ' && buf[offs]!='\t' && buf[offs]!='\n' && buf[offs]!='\r' && buf[offs]!=0 )
916            {
917                if( buf[offs]=='0' )
918                {
919                    was0 = true;
920                    offs++;
921                    continue;
922                }
923                if( buf[offs]=='1' )
924                {
925                    was1 = true;
926                    offs++;
927                    continue;
928                }
929                throw new alglib.alglibexception(emsg);
930            }
931            if( (!was0) && (!was1) )
932                throw new alglib.alglibexception(emsg);
933            if( was0 && was1 )
934                throw new alglib.alglibexception(emsg);
935            return was1 ? true : false;
936        }
937
938        /************************************************************************
939        This function serializes integer value into buffer
940
941        v           integer value to be serialized
942        buf         buffer, at least 11 characters wide
943        offs        offset in the buffer
944       
945        after return from this function, offs points to the char's past the value
946        being read.
947
948        This function raises an error in case unexpected symbol is found
949        ************************************************************************/
950        private static void int2str(int v, char[] buf, ref int offs)
951        {
952            int i;
953            byte[] _bytes = System.BitConverter.GetBytes((int)v);
954            byte[]  bytes = new byte[9];
955            int[] sixbits = new int[12];
956            byte c;
957           
958            //
959            // copy v to array of bytes, sign extending it and
960            // converting to little endian order. Additionally,
961            // we set 9th byte to zero in order to simplify
962            // conversion to six-bit representation
963            //
964            if( !System.BitConverter.IsLittleEndian )
965                System.Array.Reverse(_bytes);
966            c = v<0 ? (byte)0xFF : (byte)0x00;
967            for(i=0; i<sizeof(int); i++)
968                bytes[i] = _bytes[i];
969            for(i=sizeof(int); i<8; i++)
970                bytes[i] = c;
971            bytes[8] = 0;
972           
973            //
974            // convert to six-bit representation, output
975            //
976            // NOTE: last 12th element of sixbits is always zero, we do not output it
977            //
978            threebytes2foursixbits(bytes, 0, sixbits, 0);
979            threebytes2foursixbits(bytes, 3, sixbits, 4);
980            threebytes2foursixbits(bytes, 6, sixbits, 8);       
981            for(i=0; i<SER_ENTRY_LENGTH; i++)
982                buf[offs+i] = sixbits2char(sixbits[i]);
983            offs += SER_ENTRY_LENGTH;
984        }
985
986        /************************************************************************
987        This function unserializes integer value from string
988
989        buf         buffer which contains value; leading spaces/tabs/newlines are
990                    ignored, traling spaces/tabs/newlines are treated as  end  of
991                    the integer value.
992        offs        offset in the buffer
993       
994        after return from this function, offs points to the char's past the value
995        being read.
996
997        This function raises an error in case unexpected symbol is found
998        ************************************************************************/
999        private static int str2int(char[] buf, ref int offs)
1000        {
1001            string emsg =       "ALGLIB: unable to read integer value from stream";
1002            string emsg3264 =   "ALGLIB: unable to read integer value from stream (value does not fit into 32 bits)";
1003            int[] sixbits = new int[12];
1004            byte[] bytes = new byte[9];
1005            byte[] _bytes = new byte[sizeof(int)];
1006            int sixbitsread, i;
1007            byte c;
1008           
1009            //
1010            // 1. skip leading spaces
1011            // 2. read and decode six-bit digits
1012            // 3. set trailing digits to zeros
1013            // 4. convert to little endian 64-bit integer representation
1014            // 5. check that we fit into int
1015            // 6. convert to big endian representation, if needed
1016            //
1017            sixbitsread = 0;
1018            while( buf[offs]==' ' || buf[offs]=='\t' || buf[offs]=='\n' || buf[offs]=='\r' )
1019                offs++;
1020            while( buf[offs]!=' ' && buf[offs]!='\t' && buf[offs]!='\n' && buf[offs]!='\r' && buf[offs]!=0 )
1021            {
1022                int d;
1023                d = char2sixbits(buf[offs]);
1024                if( d<0 || sixbitsread>=SER_ENTRY_LENGTH )
1025                    throw new alglib.alglibexception(emsg);
1026                sixbits[sixbitsread] = d;
1027                sixbitsread++;
1028                offs++;
1029            }
1030            if( sixbitsread==0 )
1031                throw new alglib.alglibexception(emsg);
1032            for(i=sixbitsread; i<12; i++)
1033                sixbits[i] = 0;
1034            foursixbits2threebytes(sixbits, 0, bytes, 0);
1035            foursixbits2threebytes(sixbits, 4, bytes, 3);
1036            foursixbits2threebytes(sixbits, 8, bytes, 6);
1037            c = (bytes[sizeof(int)-1] & 0x80)!=0 ? (byte)0xFF : (byte)0x00;
1038            for(i=sizeof(int); i<8; i++)
1039                if( bytes[i]!=c )
1040                    throw new alglib.alglibexception(emsg3264);
1041            for(i=0; i<sizeof(int); i++)
1042                _bytes[i] = bytes[i];       
1043            if( !System.BitConverter.IsLittleEndian )
1044                System.Array.Reverse(_bytes);
1045            return System.BitConverter.ToInt32(_bytes,0);
1046        }   
1047       
1048       
1049        /************************************************************************
1050        This function serializes double value into buffer
1051
1052        v           double value to be serialized
1053        buf         buffer, at least 11 characters wide
1054        offs        offset in the buffer
1055       
1056        after return from this function, offs points to the char's past the value
1057        being read.
1058        ************************************************************************/
1059        private static void double2str(double v, char[] buf, ref int offs)
1060        {
1061            int i;
1062            int[] sixbits = new int[12];
1063            byte[] bytes = new byte[9];
1064
1065            //
1066            // handle special quantities
1067            //
1068            if( System.Double.IsNaN(v) )
1069            {
1070                buf[offs+0] = '.';
1071                buf[offs+1] = 'n';
1072                buf[offs+2] = 'a';
1073                buf[offs+3] = 'n';
1074                buf[offs+4] = '_';
1075                buf[offs+5] = '_';
1076                buf[offs+6] = '_';
1077                buf[offs+7] = '_';
1078                buf[offs+8] = '_';
1079                buf[offs+9] = '_';
1080                buf[offs+10] = '_';
1081                offs += SER_ENTRY_LENGTH;
1082                return;
1083            }
1084            if( System.Double.IsPositiveInfinity(v) )
1085            {
1086                buf[offs+0] = '.';
1087                buf[offs+1] = 'p';
1088                buf[offs+2] = 'o';
1089                buf[offs+3] = 's';
1090                buf[offs+4] = 'i';
1091                buf[offs+5] = 'n';
1092                buf[offs+6] = 'f';
1093                buf[offs+7] = '_';
1094                buf[offs+8] = '_';
1095                buf[offs+9] = '_';
1096                buf[offs+10] = '_';
1097                offs += SER_ENTRY_LENGTH;
1098                return;
1099            }
1100            if( System.Double.IsNegativeInfinity(v) )
1101            {
1102                buf[offs+0] = '.';
1103                buf[offs+1] = 'n';
1104                buf[offs+2] = 'e';
1105                buf[offs+3] = 'g';
1106                buf[offs+4] = 'i';
1107                buf[offs+5] = 'n';
1108                buf[offs+6] = 'f';
1109                buf[offs+7] = '_';
1110                buf[offs+8] = '_';
1111                buf[offs+9] = '_';
1112                buf[offs+10] = '_';
1113                offs += SER_ENTRY_LENGTH;
1114                return;
1115            }
1116           
1117            //
1118            // process general case:
1119            // 1. copy v to array of chars
1120            // 2. set 9th byte to zero in order to simplify conversion to six-bit representation
1121            // 3. convert to little endian (if needed)
1122            // 4. convert to six-bit representation
1123            //    (last 12th element of sixbits is always zero, we do not output it)
1124            //
1125            byte[] _bytes = System.BitConverter.GetBytes((double)v);
1126            if( !System.BitConverter.IsLittleEndian )
1127                System.Array.Reverse(_bytes);
1128            for(i=0; i<sizeof(double); i++)
1129                bytes[i] = _bytes[i];
1130            for(i=sizeof(double); i<9; i++)
1131                bytes[i] = 0;
1132            threebytes2foursixbits(bytes, 0, sixbits, 0);
1133            threebytes2foursixbits(bytes, 3, sixbits, 4);
1134            threebytes2foursixbits(bytes, 6, sixbits, 8);
1135            for(i=0; i<SER_ENTRY_LENGTH; i++)
1136                buf[offs+i] = sixbits2char(sixbits[i]);
1137            offs += SER_ENTRY_LENGTH;
1138        }
1139
1140        /************************************************************************
1141        This function unserializes double value from string
1142
1143        buf         buffer which contains value; leading spaces/tabs/newlines are
1144                    ignored, traling spaces/tabs/newlines are treated as  end  of
1145                    the double value.
1146        offs        offset in the buffer
1147       
1148        after return from this function, offs points to the char's past the value
1149        being read.
1150
1151        This function raises an error in case unexpected symbol is found
1152        ************************************************************************/
1153        private static double str2double(char[] buf, ref int offs)
1154        {
1155            string emsg = "ALGLIB: unable to read double value from stream";
1156            int[] sixbits = new int[12];
1157            byte[]  bytes = new byte[9];
1158            byte[] _bytes = new byte[sizeof(double)];
1159            int sixbitsread, i;
1160           
1161           
1162            //
1163            // skip leading spaces
1164            //
1165            while( buf[offs]==' ' || buf[offs]=='\t' || buf[offs]=='\n' || buf[offs]=='\r' )
1166                offs++;
1167           
1168             
1169            //
1170            // Handle special cases
1171            //
1172            if( buf[offs]=='.' )
1173            {
1174                string s = new string(buf, offs, SER_ENTRY_LENGTH);
1175                if( s==".nan_______" )
1176                {
1177                    offs += SER_ENTRY_LENGTH;
1178                    return System.Double.NaN;
1179                }
1180                if( s==".posinf____" )
1181                {
1182                    offs += SER_ENTRY_LENGTH;
1183                    return System.Double.PositiveInfinity;
1184                }
1185                if( s==".neginf____" )
1186                {
1187                    offs += SER_ENTRY_LENGTH;
1188                    return System.Double.NegativeInfinity;
1189                }
1190                throw new alglib.alglibexception(emsg);
1191            }
1192           
1193            //
1194            // General case:
1195            // 1. read and decode six-bit digits
1196            // 2. check that all 11 digits were read
1197            // 3. set last 12th digit to zero (needed for simplicity of conversion)
1198            // 4. convert to 8 bytes
1199            // 5. convert to big endian representation, if needed
1200            //
1201            sixbitsread = 0;
1202            while( buf[offs]!=' ' && buf[offs]!='\t' && buf[offs]!='\n' && buf[offs]!='\r' && buf[offs]!=0 )
1203            {
1204                int d;
1205                d = char2sixbits(buf[offs]);
1206                if( d<0 || sixbitsread>=SER_ENTRY_LENGTH )
1207                    throw new alglib.alglibexception(emsg);
1208                sixbits[sixbitsread] = d;
1209                sixbitsread++;
1210                offs++;
1211            }
1212            if( sixbitsread!=SER_ENTRY_LENGTH )
1213                throw new alglib.alglibexception(emsg);
1214            sixbits[SER_ENTRY_LENGTH] = 0;
1215            foursixbits2threebytes(sixbits, 0, bytes, 0);
1216            foursixbits2threebytes(sixbits, 4, bytes, 3);
1217            foursixbits2threebytes(sixbits, 8, bytes, 6);
1218            for(i=0; i<sizeof(double); i++)
1219                _bytes[i] = bytes[i];       
1220            if( !System.BitConverter.IsLittleEndian )
1221                System.Array.Reverse(_bytes);       
1222            return System.BitConverter.ToDouble(_bytes,0);
1223        }
1224    }}
Note: See TracBrowser for help on using the repository browser.