Free cookie consent management tool by TermsFeed Policy Generator

source: branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/csharp/src/ap.cs @ 15311

Last change on this file since 15311 was 15311, checked in by gkronber, 7 years ago

#2789 exploration of alglib solver for non-linear optimization with non-linear constraints

File size: 82.8 KB
Line 
1/**************************************************************************
2ALGLIB 3.11.0 (source code generated 2017-05-11)
3Copyright (c) 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    Critical failure, resilts in immediate termination of entire program.
196    ********************************************************************/
197    public static void AE_CRITICAL_ASSERT(bool x)
198    {
199        if( !x )
200            System.Environment.FailFast("ALGLIB: critical error");
201    }
202   
203    /********************************************************************
204    ALGLIB object, parent  class  for  all  internal  AlgoPascal  objects
205    managed by ALGLIB.
206   
207    Any internal AlgoPascal object inherits from this class.
208   
209    User-visible objects inherit from alglibobject (see below).
210    ********************************************************************/
211    public abstract class apobject
212    {
213        public abstract void init();
214        public abstract apobject make_copy();
215    }
216   
217    /********************************************************************
218    ALGLIB object, parent class for all user-visible objects  managed  by
219    ALGLIB.
220   
221    Methods:
222        _deallocate()       deallocation:
223                            * in managed ALGLIB it does nothing
224                            * in native ALGLIB it clears  dynamic  memory
225                              being  hold  by  object  and  sets internal
226                              reference to null.
227        make_copy()         creates deep copy of the object.
228                            Works in both managed and native versions  of
229                            ALGLIB.
230    ********************************************************************/
231    public abstract class alglibobject
232    {
233        public virtual void _deallocate() {}
234        public abstract alglibobject make_copy();
235    }
236   
237    /********************************************************************
238    Deallocation of ALGLIB object:
239    * in managed ALGLIB this method just sets refence to null
240    * in native ALGLIB call of this method:
241      1) clears dynamic memory being hold by  object  and  sets  internal
242         reference to null.
243      2) sets to null variable being passed to this method
244   
245    IMPORTANT (1): in  native  edition  of  ALGLIB,  obj becomes unusable
246                   after this call!!!  It  is  possible  to  save  a copy
247                   of reference in another variable (original variable is
248                   set to null), but any attempt to work with this object
249                   will crash your program.
250   
251    IMPORTANT (2): memory ownen by object will be recycled by GC  in  any
252                   case. This method just enforced IMMEDIATE deallocation.
253    ********************************************************************/
254    public static void deallocateimmediately<T>(ref T obj) where T : alglib.alglibobject
255    {
256        obj._deallocate();
257        obj = null;
258    }
259
260    /********************************************************************
261    Allocation counter:
262    * in managed ALGLIB it always returns 0 (dummy code)
263    * in native ALGLIB it returns current value of the allocation counter
264      (if it was activated)
265    ********************************************************************/
266    public static long alloc_counter()
267    {
268        return 0;
269    }
270   
271    /********************************************************************
272    Activization of the allocation counter:
273    * in managed ALGLIB it does nothing (dummy code)
274    * in native ALGLIB it turns on allocation counting.
275    ********************************************************************/
276    public static void alloc_counter_activate()
277    {
278    }
279   
280    /********************************************************************
281    reverse communication structure
282    ********************************************************************/
283    public class rcommstate : apobject
284    {
285        public rcommstate()
286        {
287            init();
288        }
289        public override void init()
290        {
291            stage = -1;
292            ia = new int[0];
293            ba = new bool[0];
294            ra = new double[0];
295            ca = new alglib.complex[0];
296        }
297        public override apobject make_copy()
298        {
299            rcommstate result = new rcommstate();
300            result.stage = stage;
301            result.ia = (int[])ia.Clone();
302            result.ba = (bool[])ba.Clone();
303            result.ra = (double[])ra.Clone();
304            result.ca = (alglib.complex[])ca.Clone();
305            return result;
306        }
307        public int stage;
308        public int[] ia;
309        public bool[] ba;
310        public double[] ra;
311        public alglib.complex[] ca;
312    };
313
314    /********************************************************************
315    internal functions
316    ********************************************************************/
317    public class ap
318    {
319        public static int len<T>(T[] a)
320        { return a.Length; }
321        public static int rows<T>(T[,] a)
322        { return a.GetLength(0); }
323        public static int cols<T>(T[,] a)
324        { return a.GetLength(1); }
325        public static void swap<T>(ref T a, ref T b)
326        {
327            T t = a;
328            a = b;
329            b = t;
330        }
331       
332        public static void assert(bool cond, string s)
333        {
334            if( !cond )
335                throw new alglibexception(s);
336        }
337       
338        public static void assert(bool cond)
339        {
340            assert(cond, "ALGLIB: assertion failed");
341        }
342       
343        /****************************************************************
344        returns dps (digits-of-precision) value corresponding to threshold.
345        dps(0.9)  = dps(0.5)  = dps(0.1) = 0
346        dps(0.09) = dps(0.05) = dps(0.01) = 1
347        and so on
348        ****************************************************************/
349        public static int threshold2dps(double threshold)
350        {
351            int result = 0;
352            double t;
353            for (result = 0, t = 1; t / 10 > threshold*(1+1E-10); result++, t /= 10) ;
354            return result;
355        }
356
357        /****************************************************************
358        prints formatted complex
359        ****************************************************************/
360        public static string format(complex a, int _dps)
361        {
362            int dps = Math.Abs(_dps);
363            string fmt = _dps>=0 ? "F" : "E";
364            string fmtx = String.Format("{{0:"+fmt+"{0}}}", dps);
365            string fmty = String.Format("{{0:"+fmt+"{0}}}", dps);
366            string result = String.Format(fmtx, a.x) + (a.y >= 0 ? "+" : "-") + String.Format(fmty, Math.Abs(a.y)) + "i";
367            result = result.Replace(',', '.');
368            return result;
369        }
370
371        /****************************************************************
372        prints formatted array
373        ****************************************************************/
374        public static string format(bool[] a)
375        {
376            string[] result = new string[len(a)];
377            int i;
378            for(i=0; i<len(a); i++)
379                if( a[i] )
380                    result[i] = "true";
381                else
382                    result[i] = "false";
383            return "{"+String.Join(",",result)+"}";
384        }
385       
386        /****************************************************************
387        prints formatted array
388        ****************************************************************/
389        public static string format(int[] a)
390        {
391            string[] result = new string[len(a)];
392            int i;
393            for (i = 0; i < len(a); i++)
394                result[i] = a[i].ToString();
395            return "{" + String.Join(",", result) + "}";
396        }
397
398        /****************************************************************
399        prints formatted array
400        ****************************************************************/
401        public static string format(double[] a, int _dps)
402        {
403            int dps = Math.Abs(_dps);
404            string sfmt = _dps >= 0 ? "F" : "E";
405            string fmt = String.Format("{{0:" + sfmt + "{0}}}", dps);
406            string[] result = new string[len(a)];
407            int i;
408            for (i = 0; i < len(a); i++)
409            {
410                result[i] = String.Format(fmt, a[i]);
411                result[i] = result[i].Replace(',', '.');
412            }
413            return "{" + String.Join(",", result) + "}";
414        }
415
416        /****************************************************************
417        prints formatted array
418        ****************************************************************/
419        public static string format(complex[] a, int _dps)
420        {
421            int dps = Math.Abs(_dps);
422            string fmt = _dps >= 0 ? "F" : "E";
423            string fmtx = String.Format("{{0:"+fmt+"{0}}}", dps);
424            string fmty = String.Format("{{0:"+fmt+"{0}}}", dps);
425            string[] result = new string[len(a)];
426            int i;
427            for (i = 0; i < len(a); i++)
428            {
429                result[i] = String.Format(fmtx, a[i].x) + (a[i].y >= 0 ? "+" : "-") + String.Format(fmty, Math.Abs(a[i].y)) + "i";
430                result[i] = result[i].Replace(',', '.');
431            }
432            return "{" + String.Join(",", result) + "}";
433        }
434
435        /****************************************************************
436        prints formatted matrix
437        ****************************************************************/
438        public static string format(bool[,] a)
439        {
440            int i, j, m, n;
441            n = cols(a);
442            m = rows(a);
443            bool[] line = new bool[n];
444            string[] result = new string[m];
445            for (i = 0; i < m; i++)
446            {
447                for (j = 0; j < n; j++)
448                    line[j] = a[i, j];
449                result[i] = format(line);
450            }
451            return "{" + String.Join(",", result) + "}";
452        }
453
454        /****************************************************************
455        prints formatted matrix
456        ****************************************************************/
457        public static string format(int[,] a)
458        {
459            int i, j, m, n;
460            n = cols(a);
461            m = rows(a);
462            int[] line = new int[n];
463            string[] result = new string[m];
464            for (i = 0; i < m; i++)
465            {
466                for (j = 0; j < n; j++)
467                    line[j] = a[i, j];
468                result[i] = format(line);
469            }
470            return "{" + String.Join(",", result) + "}";
471        }
472
473        /****************************************************************
474        prints formatted matrix
475        ****************************************************************/
476        public static string format(double[,] a, int dps)
477        {
478            int i, j, m, n;
479            n = cols(a);
480            m = rows(a);
481            double[] line = new double[n];
482            string[] result = new string[m];
483            for (i = 0; i < m; i++)
484            {
485                for (j = 0; j < n; j++)
486                    line[j] = a[i, j];
487                result[i] = format(line, dps);
488            }
489            return "{" + String.Join(",", result) + "}";
490        }
491
492        /****************************************************************
493        prints formatted matrix
494        ****************************************************************/
495        public static string format(complex[,] a, int dps)
496        {
497            int i, j, m, n;
498            n = cols(a);
499            m = rows(a);
500            complex[] line = new complex[n];
501            string[] result = new string[m];
502            for (i = 0; i < m; i++)
503            {
504                for (j = 0; j < n; j++)
505                    line[j] = a[i, j];
506                result[i] = format(line, dps);
507            }
508            return "{" + String.Join(",", result) + "}";
509        }
510
511        /****************************************************************
512        checks that matrix is symmetric.
513        max|A-A^T| is calculated; if it is within 1.0E-14 of max|A|,
514        matrix is considered symmetric
515        ****************************************************************/
516        public static bool issymmetric(double[,] a)
517        {
518            int i, j, n;
519            double err, mx, v1, v2;
520            if( rows(a)!=cols(a) )
521                return false;
522            n = rows(a);
523            if( n==0 )
524                return true;
525            mx = 0;
526            err = 0;
527            for( i=0; i<n; i++)
528            {
529                for(j=i+1; j<n; j++)
530                {
531                    v1 = a[i,j];
532                    v2 = a[j,i];
533                    if( !math.isfinite(v1) )
534                        return false;
535                    if( !math.isfinite(v2) )
536                        return false;
537                    err = Math.Max(err, Math.Abs(v1-v2));
538                    mx  = Math.Max(mx,  Math.Abs(v1));
539                    mx  = Math.Max(mx,  Math.Abs(v2));
540                }
541                v1 = a[i,i];
542                if( !math.isfinite(v1) )
543                    return false;
544                mx = Math.Max(mx, Math.Abs(v1));
545            }
546            if( mx==0 )
547                return true;
548            return err/mx<=1.0E-14;
549        }
550       
551        /****************************************************************
552        checks that matrix is Hermitian.
553        max|A-A^H| is calculated; if it is within 1.0E-14 of max|A|,
554        matrix is considered Hermitian
555        ****************************************************************/
556        public static bool ishermitian(complex[,] a)
557        {
558            int i, j, n;
559            double err, mx;
560            complex v1, v2, vt;
561            if( rows(a)!=cols(a) )
562                return false;
563            n = rows(a);
564            if( n==0 )
565                return true;
566            mx = 0;
567            err = 0;
568            for( i=0; i<n; i++)
569            {
570                for(j=i+1; j<n; j++)
571                {
572                    v1 = a[i,j];
573                    v2 = a[j,i];
574                    if( !math.isfinite(v1.x) )
575                        return false;
576                    if( !math.isfinite(v1.y) )
577                        return false;
578                    if( !math.isfinite(v2.x) )
579                        return false;
580                    if( !math.isfinite(v2.y) )
581                        return false;
582                    vt.x = v1.x-v2.x;
583                    vt.y = v1.y+v2.y;
584                    err = Math.Max(err, math.abscomplex(vt));
585                    mx  = Math.Max(mx,  math.abscomplex(v1));
586                    mx  = Math.Max(mx,  math.abscomplex(v2));
587                }
588                v1 = a[i,i];
589                if( !math.isfinite(v1.x) )
590                    return false;
591                if( !math.isfinite(v1.y) )
592                    return false;
593                err = Math.Max(err, Math.Abs(v1.y));
594                mx = Math.Max(mx, math.abscomplex(v1));
595            }
596            if( mx==0 )
597                return true;
598            return err/mx<=1.0E-14;
599        }
600       
601       
602        /****************************************************************
603        Forces symmetricity by copying upper half of A to the lower one
604        ****************************************************************/
605        public static bool forcesymmetric(double[,] a)
606        {
607            int i, j, n;
608            if( rows(a)!=cols(a) )
609                return false;
610            n = rows(a);
611            if( n==0 )
612                return true;
613            for( i=0; i<n; i++)
614                for(j=i+1; j<n; j++)
615                    a[i,j] = a[j,i];
616            return true;
617        }
618       
619        /****************************************************************
620        Forces Hermiticity by copying upper half of A to the lower one
621        ****************************************************************/
622        public static bool forcehermitian(complex[,] a)
623        {
624            int i, j, n;
625            complex v;
626            if( rows(a)!=cols(a) )
627                return false;
628            n = rows(a);
629            if( n==0 )
630                return true;
631            for( i=0; i<n; i++)
632                for(j=i+1; j<n; j++)
633                {
634                    v = a[j,i];
635                    a[i,j].x = v.x;
636                    a[i,j].y = -v.y;
637                }
638            return true;
639        }
640    };
641   
642    /********************************************************************
643    math functions
644    ********************************************************************/
645    public class math
646    {
647        //public static System.Random RndObject = new System.Random(System.DateTime.Now.Millisecond);
648        public static System.Random rndobject = new System.Random(System.DateTime.Now.Millisecond + 1000*System.DateTime.Now.Second + 60*1000*System.DateTime.Now.Minute);
649
650        public const double machineepsilon = 5E-16;
651        public const double maxrealnumber = 1E300;
652        public const double minrealnumber = 1E-300;
653       
654        public static bool isfinite(double d)
655        {
656            return !System.Double.IsNaN(d) && !System.Double.IsInfinity(d);
657        }
658       
659        public static double randomreal()
660        {
661            double r = 0;
662            lock(rndobject){ r = rndobject.NextDouble(); }
663            return r;
664        }
665        public static int randominteger(int N)
666        {
667            int r = 0;
668            lock(rndobject){ r = rndobject.Next(N); }
669            return r;
670        }
671        public static double sqr(double X)
672        {
673            return X*X;
674        }       
675        public static double abscomplex(complex z)
676        {
677            double w;
678            double xabs;
679            double yabs;
680            double v;
681   
682            xabs = System.Math.Abs(z.x);
683            yabs = System.Math.Abs(z.y);
684            w = xabs>yabs ? xabs : yabs;
685            v = xabs<yabs ? xabs : yabs;
686            if( v==0 )
687                return w;
688            else
689            {
690                double t = v/w;
691                return w*System.Math.Sqrt(1+t*t);
692            }
693        }
694        public static complex conj(complex z)
695        {
696            return new complex(z.x, -z.y);
697        }   
698        public static complex csqr(complex z)
699        {
700            return new complex(z.x*z.x-z.y*z.y, 2*z.x*z.y);
701        }
702
703    }
704
705    /*
706     * CSV functionality
707     */
708     
709    public static int CSV_DEFAULT      = 0x0;
710    public static int CSV_SKIP_HEADERS = 0x1;
711   
712    /*
713     * CSV operations: reading CSV file to real matrix.
714     *
715     * This function reads CSV  file  and  stores  its  contents  to  double
716     * precision 2D array. Format of the data file must conform to RFC  4180
717     * specification, with additional notes:
718     * - file size should be less than 2GB
719     * - ASCI encoding, UTF-8 without BOM (in header names) are supported
720     * - any character (comma/tab/space) may be used as field separator,  as
721     *   long as it is distinct from one used for decimal point
722     * - multiple subsequent field separators (say, two  spaces) are treated
723     *   as MULTIPLE separators, not one big separator
724     * - both comma and full stop may be used as decimal point. Parser  will
725     *   automatically determine specific character being used.  Both  fixed
726     *   and exponential number formats are  allowed.   Thousand  separators
727     *   are NOT allowed.
728     * - line may end with \n (Unix style) or \r\n (Windows  style),  parser
729     *   will automatically adapt to chosen convention
730     * - escaped fields (ones in double quotes) are not supported
731     *
732     * INPUT PARAMETERS:
733     *     filename        relative/absolute path
734     *     separator       character used to separate fields.  May  be  ' ',
735     *                     ',', '\t'. Other separators are possible too.
736     *     flags           several values combined with bitwise OR:
737     *                     * alglib::CSV_SKIP_HEADERS -  if present, first row
738     *                       contains headers  and  will  be  skipped.   Its
739     *                       contents is used to determine fields count, and
740     *                       that's all.
741     *                     If no flags are specified, default value 0x0  (or
742     *                     alglib::CSV_DEFAULT, which is same) should be used.
743     *                     
744     * OUTPUT PARAMETERS:
745     *     out             2D matrix, CSV file parsed with atof()
746     *     
747     * HANDLING OF SPECIAL CASES:
748     * - file does not exist - alglib::ap_error exception is thrown
749     * - empty file - empty array is returned (no exception)
750     * - skip_first_row=true, only one row in file - empty array is returned
751     * - field contents is not recognized by atof() - field value is replaced
752     *   by 0.0
753     */
754    public static void read_csv(string filename, char separator, int flags, out double[,] matrix)
755    {
756        //
757        // Parameters
758        //
759        bool skip_first_row = (flags&CSV_SKIP_HEADERS)!=0;
760       
761        //
762        // Prepare empty output array
763        //
764        matrix = new double[0,0];
765       
766        //
767        // Read file, normalize file contents:
768        // * replace 0x0 by spaces
769        // * remove trailing spaces and newlines
770        // * append trailing '\n' and '\0' characters
771        // Return if file contains only spaces/newlines.
772        //
773        byte b_space = System.Convert.ToByte(' ');
774        byte b_tab   = System.Convert.ToByte('\t');
775        byte b_lf    = System.Convert.ToByte('\n');
776        byte b_cr    = System.Convert.ToByte('\r');
777        byte b_comma = System.Convert.ToByte(',');
778        byte b_fullstop= System.Convert.ToByte('.');
779        byte[] v0 = System.IO.File.ReadAllBytes(filename);
780        if( v0.Length==0 )
781            return;
782        byte[] v1 = new byte[v0.Length+2];
783        int filesize = v0.Length;
784        for(int i=0; i<filesize; i++)
785            v1[i] = v0[i]==0 ? b_space : v0[i];
786        for(; filesize>0; )
787        {
788            byte c = v1[filesize-1];
789            if( c==b_space || c==b_tab || c==b_cr || c==b_lf )
790            {
791                filesize--;
792                continue;
793            }
794            break;
795        }
796        if( filesize==0 )
797            return;
798        v1[filesize+0] = b_lf;
799        v1[filesize+1] = 0x0;
800        filesize+=2;
801       
802       
803        //
804        // Scan dataset.
805        //
806        int rows_count, cols_count, max_length = 0;
807        cols_count = 1;
808        for(int idx=0; idx<filesize; idx++)
809        {
810            if( v1[idx]==separator )
811                cols_count++;
812            if( v1[idx]==b_lf )
813                break;
814        }
815        rows_count = 0;
816        for(int idx=0; idx<filesize; idx++)
817            if( v1[idx]==b_lf )
818                rows_count++;
819        if( rows_count==1 && skip_first_row ) // empty output, return
820            return;
821        int[] offsets = new int[rows_count*cols_count];
822        int[] lengths = new int[rows_count*cols_count];
823        int cur_row_idx = 0;
824        for(int row_start=0; v1[row_start]!=0x0; )
825        {
826            // determine row length
827            int row_length;
828            for(row_length=0; v1[row_start+row_length]!=b_lf; row_length++);
829           
830            // determine cols count, perform integrity check
831            int cur_cols_cnt=1;
832            for(int idx=0; idx<row_length; idx++)
833                if( v1[row_start+idx]==separator )
834                    cur_cols_cnt++;
835            if( cols_count!=cur_cols_cnt )
836                throw new alglib.alglibexception("read_csv: non-rectangular contents, rows have different sizes");
837           
838            // store offsets and lengths of the fields
839            int cur_offs = 0;
840            int cur_col_idx = 0;
841            for(int idx=0; idx<row_length+1; idx++)
842                if( v1[row_start+idx]==separator || v1[row_start+idx]==b_lf )
843                {
844                    offsets[cur_row_idx*cols_count+cur_col_idx] = row_start+cur_offs;
845                    lengths[cur_row_idx*cols_count+cur_col_idx] = idx-cur_offs;
846                    max_length = idx-cur_offs>max_length ? idx-cur_offs : max_length;
847                    cur_offs = idx+1;
848                    cur_col_idx++;
849                }
850           
851            // advance row start
852            cur_row_idx++;
853            row_start = row_start+row_length+1;
854        }
855       
856        //
857        // Convert
858        //
859        int row0 = skip_first_row ? 1 : 0;
860        int row1 = rows_count;
861        System.Globalization.CultureInfo culture = System.Globalization.CultureInfo.CreateSpecificCulture(""); // invariant culture
862        matrix = new double[row1-row0, cols_count];
863        alglib.AE_CRITICAL_ASSERT(culture.NumberFormat.NumberDecimalSeparator==".");
864        for(int ridx=row0; ridx<row1; ridx++)
865            for(int cidx=0; cidx<cols_count; cidx++)
866            {
867                int field_len  = lengths[ridx*cols_count+cidx];
868                int field_offs = offsets[ridx*cols_count+cidx];
869               
870                // replace , by full stop
871                for(int idx=0; idx<field_len; idx++)
872                    if( v1[field_offs+idx]==b_comma )
873                        v1[field_offs+idx] = b_fullstop;
874               
875                // convert
876                string s_val = System.Text.Encoding.ASCII.GetString(v1, field_offs, field_len);
877                double d_val;
878                Double.TryParse(s_val, System.Globalization.NumberStyles.Float, culture, out d_val);
879                matrix[ridx-row0,cidx] = d_val;
880            }
881    }
882   
883   
884    /********************************************************************
885    serializer object (should not be used directly)
886    ********************************************************************/
887    public class serializer
888    {
889        enum SMODE { DEFAULT, ALLOC, TO_STRING, FROM_STRING, TO_STREAM, FROM_STREAM };
890        private const int SER_ENTRIES_PER_ROW = 5;
891        private const int SER_ENTRY_LENGTH    = 11;
892       
893        private SMODE mode;
894        private int entries_needed;
895        private int entries_saved;
896        private int bytes_asked;
897        private int bytes_written;
898        private int bytes_read;
899        private char[] out_str;
900        private char[] in_str;
901        private System.IO.Stream io_stream;
902       
903        // local temporaries
904        private char[] entry_buf_char;
905        private byte[] entry_buf_byte;
906       
907        public serializer()
908        {
909            mode = SMODE.DEFAULT;
910            entries_needed = 0;
911            bytes_asked = 0;
912            entry_buf_byte = new byte[SER_ENTRY_LENGTH+2];
913            entry_buf_char = new char[SER_ENTRY_LENGTH+2];
914        }
915
916        public void clear_buffers()
917        {
918            out_str = null;
919            in_str = null;
920            io_stream = null;
921        }
922
923        public void alloc_start()
924        {
925            entries_needed = 0;
926            bytes_asked = 0;
927            mode = SMODE.ALLOC;
928        }
929
930        public void alloc_entry()
931        {
932            if( mode!=SMODE.ALLOC )
933                throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");
934            entries_needed++;
935        }
936
937        private int get_alloc_size()
938        {
939            int rows, lastrowsize, result;
940           
941            // check and change mode
942            if( mode!=SMODE.ALLOC )
943                throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");
944           
945            // if no entries needes (degenerate case)
946            if( entries_needed==0 )
947            {
948                bytes_asked = 4; /* a pair of chars for \r\n, one for space, one for dot */
949                return bytes_asked;
950            }
951           
952            // non-degenerate case
953            rows = entries_needed/SER_ENTRIES_PER_ROW;
954            lastrowsize = SER_ENTRIES_PER_ROW;
955            if( entries_needed%SER_ENTRIES_PER_ROW!=0 )
956            {
957                lastrowsize = entries_needed%SER_ENTRIES_PER_ROW;
958                rows++;
959            }
960           
961            // calculate result size
962            result  = ((rows-1)*SER_ENTRIES_PER_ROW+lastrowsize)*SER_ENTRY_LENGTH;  /* data size */
963            result +=  (rows-1)*(SER_ENTRIES_PER_ROW-1)+(lastrowsize-1);            /* space symbols */
964            result += rows*2;                                                       /* newline symbols */
965            result += 1;                                                            /* trailing dot */
966            bytes_asked = result;
967            return result;
968        }
969
970        public void sstart_str()
971        {
972            int allocsize = get_alloc_size();
973           
974            // clear input/output buffers which may hold pointers to unneeded memory
975            // NOTE: it also helps us to avoid errors when data are written to incorrect location
976            clear_buffers();
977           
978            // check and change mode
979            if( mode!=SMODE.ALLOC )
980                throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");
981            mode = SMODE.TO_STRING;
982           
983            // other preparations
984            out_str = new char[allocsize];
985            entries_saved = 0;
986            bytes_written = 0;
987        }
988
989        public void sstart_stream(System.IO.Stream o_stream)
990        {   
991            // clear input/output buffers which may hold pointers to unneeded memory
992            // NOTE: it also helps us to avoid errors when data are written to incorrect location
993            clear_buffers();
994           
995            // check and change mode
996            if( mode!=SMODE.ALLOC )
997                throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");
998            mode = SMODE.TO_STREAM;
999            io_stream = o_stream;
1000        }
1001
1002        public void ustart_str(string s)
1003        {
1004            // clear input/output buffers which may hold pointers to unneeded memory
1005            // NOTE: it also helps us to avoid errors when data are written to incorrect location
1006            clear_buffers();
1007           
1008            // check and change mode
1009            if( mode!=SMODE.DEFAULT )
1010                throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");
1011            mode = SMODE.FROM_STRING;
1012           
1013            in_str = s.ToCharArray();
1014            bytes_read = 0;
1015        }
1016
1017        public void ustart_stream(System.IO.Stream i_stream)
1018        {
1019            // clear input/output buffers which may hold pointers to unneeded memory
1020            // NOTE: it also helps us to avoid errors when data are written to incorrect location
1021            clear_buffers();
1022           
1023            // check and change mode
1024            if( mode!=SMODE.DEFAULT )
1025                throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");
1026            mode = SMODE.FROM_STREAM;
1027            io_stream = i_stream;
1028        }
1029
1030        private void serialize_value(bool v0, int v1, double v2, int val_idx)
1031        {
1032            // prepare serialization
1033            char[] arr_out = null;
1034            int cnt_out = 0;
1035            if( mode==SMODE.TO_STRING )
1036            {
1037                arr_out = out_str;
1038                cnt_out = bytes_written;
1039            }
1040            else if( mode==SMODE.TO_STREAM )
1041            {
1042                arr_out = entry_buf_char;
1043                cnt_out = 0;
1044            }
1045            else
1046                throw new alglib.alglibexception("ALGLIB: internal error during serialization");
1047           
1048            // serialize
1049            if( val_idx==0 )
1050                bool2str(  v0, arr_out, ref cnt_out);
1051            else if( val_idx==1 )
1052                int2str(   v1, arr_out, ref cnt_out);
1053            else if( val_idx==2 )
1054                double2str(v2, arr_out, ref cnt_out);
1055            else
1056                throw new alglib.alglibexception("ALGLIB: internal error during serialization");
1057            entries_saved++;
1058            if( entries_saved%SER_ENTRIES_PER_ROW!=0 )
1059            {
1060                arr_out[cnt_out] = ' ';
1061                cnt_out++;
1062            }
1063            else
1064            {
1065                arr_out[cnt_out+0] = '\r';
1066                arr_out[cnt_out+1] = '\n';
1067                cnt_out+=2;
1068            }
1069           
1070            // post-process
1071            if( mode==SMODE.TO_STRING )
1072            {
1073                bytes_written = cnt_out;
1074                return;
1075            }
1076            else if( mode==SMODE.TO_STREAM )
1077            {
1078                for(int k=0; k<cnt_out; k++)
1079                    entry_buf_byte[k] = (byte)entry_buf_char[k];
1080                io_stream.Write(entry_buf_byte, 0, cnt_out);
1081                return;
1082            }
1083            else
1084                throw new alglib.alglibexception("ALGLIB: internal error during serialization");
1085        }
1086
1087        private void unstream_entry_char()
1088        {
1089            if( mode!=SMODE.FROM_STREAM )
1090                throw new alglib.alglibexception("ALGLIB: internal error during unserialization");
1091            int c;
1092            for(;;)
1093            {
1094                c = io_stream.ReadByte();
1095                if( c<0 )
1096                    throw new alglib.alglibexception("ALGLIB: internal error during unserialization");
1097                if( c!=' ' && c!='\t' && c!='\n' && c!='\r' )
1098                    break;
1099            }
1100            entry_buf_char[0] = (char)c;
1101            for(int k=1; k<SER_ENTRY_LENGTH; k++)
1102            {
1103                c = io_stream.ReadByte();
1104                entry_buf_char[k] = (char)c;
1105                if( c<0 || c==' ' || c=='\t' || c=='\n' || c=='\r' )
1106                    throw new alglib.alglibexception("ALGLIB: internal error during unserialization");
1107            }
1108            entry_buf_char[SER_ENTRY_LENGTH] = (char)0;
1109        }
1110
1111        public void serialize_bool(bool v)
1112        {
1113            serialize_value(v, 0, 0, 0);
1114        }
1115
1116        public void serialize_int(int v)
1117        {
1118            serialize_value(false, v, 0, 1);
1119        }
1120
1121        public void serialize_double(double v)
1122        {
1123            serialize_value(false, 0, v, 2);
1124        }
1125
1126        public bool unserialize_bool()
1127        {
1128            if( mode==SMODE.FROM_STRING )
1129                return str2bool(in_str, ref bytes_read);
1130            if( mode==SMODE.FROM_STREAM )
1131            {
1132                unstream_entry_char();
1133                int dummy = 0;
1134                return str2bool(entry_buf_char, ref dummy);
1135            }
1136            throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");
1137        }
1138
1139        public int unserialize_int()
1140        {
1141            if( mode==SMODE.FROM_STRING )
1142                return str2int(in_str, ref bytes_read);
1143            if( mode==SMODE.FROM_STREAM )
1144            {
1145                unstream_entry_char();
1146                int dummy = 0;
1147                return str2int(entry_buf_char, ref dummy);
1148            }
1149            throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");
1150        }
1151
1152        public double unserialize_double()
1153        {
1154            if( mode==SMODE.FROM_STRING )
1155                return str2double(in_str, ref bytes_read);
1156            if( mode==SMODE.FROM_STREAM )
1157            {
1158                unstream_entry_char();
1159                int dummy = 0;
1160                return str2double(entry_buf_char, ref dummy);
1161            }
1162            throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");
1163        }
1164
1165        public void stop()
1166        {
1167            if( mode==SMODE.TO_STRING )
1168            {
1169                out_str[bytes_written] = '.';
1170                bytes_written++;
1171                return;
1172            }
1173            if( mode==SMODE.FROM_STRING )
1174            {
1175                //
1176                // because input string may be from pre-3.11 serializer,
1177                // which does not include trailing dot, we do not test
1178                // string for presence of "." symbol. Anyway, because string
1179                // is not stream, we do not have to read ALL trailing symbols.
1180                //
1181                return;
1182            }
1183            if( mode==SMODE.TO_STREAM )
1184            {
1185                io_stream.WriteByte((byte)'.');
1186                return;
1187            }
1188            if( mode==SMODE.FROM_STREAM )
1189            {
1190                for(;;)
1191                {
1192                    int c = io_stream.ReadByte();
1193                    if( c==' ' || c=='\t' || c=='\n' || c=='\r' )
1194                        continue;
1195                    if( c=='.' )
1196                        break;
1197                    throw new alglib.alglibexception("ALGLIB: internal error during unserialization");
1198                }
1199                return;
1200            }
1201            throw new alglib.alglibexception("ALGLIB: internal error during unserialization");
1202        }
1203
1204        public string get_string()
1205        {
1206            if( mode!=SMODE.TO_STRING )
1207                throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");
1208             return new string(out_str, 0, bytes_written);
1209        }
1210
1211
1212        /************************************************************************
1213        This function converts six-bit value (from 0 to 63)  to  character  (only
1214        digits, lowercase and uppercase letters, minus and underscore are used).
1215
1216        If v is negative or greater than 63, this function returns '?'.
1217        ************************************************************************/
1218        private static char[] _sixbits2char_tbl = new char[64]{
1219                '0', '1', '2', '3', '4', '5', '6', '7',
1220                '8', '9', 'A', 'B', 'C', 'D', 'E', 'F',
1221                'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N',
1222                'O', 'P', 'Q', 'R', 'S', 'T', 'U', 'V',
1223                'W', 'X', 'Y', 'Z', 'a', 'b', 'c', 'd',
1224                'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l',
1225                'm', 'n', 'o', 'p', 'q', 'r', 's', 't',
1226                'u', 'v', 'w', 'x', 'y', 'z', '-', '_' };
1227        private static char sixbits2char(int v)
1228        {
1229            if( v<0 || v>63 )
1230                return '?';
1231            return _sixbits2char_tbl[v];
1232        }
1233       
1234        /************************************************************************
1235        This function converts character to six-bit value (from 0 to 63).
1236
1237        This function is inverse of ae_sixbits2char()
1238        If c is not correct character, this function returns -1.
1239        ************************************************************************/
1240        private static int[] _char2sixbits_tbl = new int[128] {
1241            -1, -1, -1, -1, -1, -1, -1, -1,
1242            -1, -1, -1, -1, -1, -1, -1, -1,
1243            -1, -1, -1, -1, -1, -1, -1, -1,
1244            -1, -1, -1, -1, -1, -1, -1, -1,
1245            -1, -1, -1, -1, -1, -1, -1, -1,
1246            -1, -1, -1, -1, -1, 62, -1, -1,
1247             0,  1,  2,  3,  4,  5,  6,  7,
1248             8,  9, -1, -1, -1, -1, -1, -1,
1249            -1, 10, 11, 12, 13, 14, 15, 16,
1250            17, 18, 19, 20, 21, 22, 23, 24,
1251            25, 26, 27, 28, 29, 30, 31, 32,
1252            33, 34, 35, -1, -1, -1, -1, 63,
1253            -1, 36, 37, 38, 39, 40, 41, 42,
1254            43, 44, 45, 46, 47, 48, 49, 50,
1255            51, 52, 53, 54, 55, 56, 57, 58,
1256            59, 60, 61, -1, -1, -1, -1, -1 };
1257        private static int char2sixbits(char c)
1258        {
1259            return (c>=0 && c<127) ? _char2sixbits_tbl[c] : -1;
1260        }
1261       
1262        /************************************************************************
1263        This function converts three bytes (24 bits) to four six-bit values
1264        (24 bits again).
1265
1266        src         array
1267        src_offs    offset of three-bytes chunk
1268        dst         array for ints
1269        dst_offs    offset of four-ints chunk
1270        ************************************************************************/
1271        private static void threebytes2foursixbits(byte[] src, int src_offs, int[] dst, int dst_offs)
1272        {
1273            dst[dst_offs+0] =  src[src_offs+0] & 0x3F;
1274            dst[dst_offs+1] = (src[src_offs+0]>>6) | ((src[src_offs+1]&0x0F)<<2);
1275            dst[dst_offs+2] = (src[src_offs+1]>>4) | ((src[src_offs+2]&0x03)<<4);
1276            dst[dst_offs+3] =  src[src_offs+2]>>2;
1277        }
1278
1279        /************************************************************************
1280        This function converts four six-bit values (24 bits) to three bytes
1281        (24 bits again).
1282
1283        src         pointer to four ints
1284        src_offs    offset of the chunk
1285        dst         pointer to three bytes
1286        dst_offs    offset of the chunk
1287        ************************************************************************/
1288        private static void foursixbits2threebytes(int[] src, int src_offs, byte[] dst, int dst_offs)
1289        {
1290            dst[dst_offs+0] =      (byte)(src[src_offs+0] | ((src[src_offs+1]&0x03)<<6));
1291            dst[dst_offs+1] = (byte)((src[src_offs+1]>>2) | ((src[src_offs+2]&0x0F)<<4));
1292            dst[dst_offs+2] = (byte)((src[src_offs+2]>>4) |  (src[src_offs+3]<<2));
1293        }
1294
1295        /************************************************************************
1296        This function serializes boolean value into buffer
1297
1298        v           boolean value to be serialized
1299        buf         buffer, at least 11 characters wide
1300        offs        offset in the buffer
1301       
1302        after return from this function, offs points to the char's past the value
1303        being read.
1304        ************************************************************************/
1305        private static void bool2str(bool v, char[] buf, ref int offs)
1306        {
1307            char c = v ? '1' : '0';
1308            int i;
1309            for(i=0; i<SER_ENTRY_LENGTH; i++)
1310                buf[offs+i] = c;
1311            offs += SER_ENTRY_LENGTH;
1312        }
1313
1314        /************************************************************************
1315        This function unserializes boolean value from buffer
1316
1317        buf         buffer which contains value; leading spaces/tabs/newlines are
1318                    ignored, traling spaces/tabs/newlines are treated as  end  of
1319                    the boolean value.
1320        offs        offset in the buffer
1321       
1322        after return from this function, offs points to the char's past the value
1323        being read.
1324
1325        This function raises an error in case unexpected symbol is found
1326        ************************************************************************/
1327        private static bool str2bool(char[] buf, ref int offs)
1328        {
1329            bool was0, was1;
1330            string emsg = "ALGLIB: unable to read boolean value from stream";
1331           
1332            was0 = false;
1333            was1 = false;
1334            while( buf[offs]==' ' || buf[offs]=='\t' || buf[offs]=='\n' || buf[offs]=='\r' )
1335                offs++;
1336            while( buf[offs]!=' ' && buf[offs]!='\t' && buf[offs]!='\n' && buf[offs]!='\r' && buf[offs]!=0 )
1337            {
1338                if( buf[offs]=='0' )
1339                {
1340                    was0 = true;
1341                    offs++;
1342                    continue;
1343                }
1344                if( buf[offs]=='1' )
1345                {
1346                    was1 = true;
1347                    offs++;
1348                    continue;
1349                }
1350                throw new alglib.alglibexception(emsg);
1351            }
1352            if( (!was0) && (!was1) )
1353                throw new alglib.alglibexception(emsg);
1354            if( was0 && was1 )
1355                throw new alglib.alglibexception(emsg);
1356            return was1 ? true : false;
1357        }
1358
1359        /************************************************************************
1360        This function serializes integer value into buffer
1361
1362        v           integer value to be serialized
1363        buf         buffer, at least 11 characters wide
1364        offs        offset in the buffer
1365       
1366        after return from this function, offs points to the char's past the value
1367        being read.
1368
1369        This function raises an error in case unexpected symbol is found
1370        ************************************************************************/
1371        private static void int2str(int v, char[] buf, ref int offs)
1372        {
1373            int i;
1374            byte[] _bytes = System.BitConverter.GetBytes((int)v);
1375            byte[]  bytes = new byte[9];
1376            int[] sixbits = new int[12];
1377            byte c;
1378           
1379            //
1380            // copy v to array of bytes, sign extending it and
1381            // converting to little endian order. Additionally,
1382            // we set 9th byte to zero in order to simplify
1383            // conversion to six-bit representation
1384            //
1385            if( !System.BitConverter.IsLittleEndian )
1386                System.Array.Reverse(_bytes);
1387            c = v<0 ? (byte)0xFF : (byte)0x00;
1388            for(i=0; i<sizeof(int); i++)
1389                bytes[i] = _bytes[i];
1390            for(i=sizeof(int); i<8; i++)
1391                bytes[i] = c;
1392            bytes[8] = 0;
1393           
1394            //
1395            // convert to six-bit representation, output
1396            //
1397            // NOTE: last 12th element of sixbits is always zero, we do not output it
1398            //
1399            threebytes2foursixbits(bytes, 0, sixbits, 0);
1400            threebytes2foursixbits(bytes, 3, sixbits, 4);
1401            threebytes2foursixbits(bytes, 6, sixbits, 8);       
1402            for(i=0; i<SER_ENTRY_LENGTH; i++)
1403                buf[offs+i] = sixbits2char(sixbits[i]);
1404            offs += SER_ENTRY_LENGTH;
1405        }
1406
1407        /************************************************************************
1408        This function unserializes integer value from string
1409
1410        buf         buffer which contains value; leading spaces/tabs/newlines are
1411                    ignored, traling spaces/tabs/newlines are treated as  end  of
1412                    the integer value.
1413        offs        offset in the buffer
1414       
1415        after return from this function, offs points to the char's past the value
1416        being read.
1417
1418        This function raises an error in case unexpected symbol is found
1419        ************************************************************************/
1420        private static int str2int(char[] buf, ref int offs)
1421        {
1422            string emsg =       "ALGLIB: unable to read integer value from stream";
1423            string emsg3264 =   "ALGLIB: unable to read integer value from stream (value does not fit into 32 bits)";
1424            int[] sixbits = new int[12];
1425            byte[] bytes = new byte[9];
1426            byte[] _bytes = new byte[sizeof(int)];
1427            int sixbitsread, i;
1428            byte c;
1429           
1430            //
1431            // 1. skip leading spaces
1432            // 2. read and decode six-bit digits
1433            // 3. set trailing digits to zeros
1434            // 4. convert to little endian 64-bit integer representation
1435            // 5. check that we fit into int
1436            // 6. convert to big endian representation, if needed
1437            //
1438            sixbitsread = 0;
1439            while( buf[offs]==' ' || buf[offs]=='\t' || buf[offs]=='\n' || buf[offs]=='\r' )
1440                offs++;
1441            while( buf[offs]!=' ' && buf[offs]!='\t' && buf[offs]!='\n' && buf[offs]!='\r' && buf[offs]!=0 )
1442            {
1443                int d;
1444                d = char2sixbits(buf[offs]);
1445                if( d<0 || sixbitsread>=SER_ENTRY_LENGTH )
1446                    throw new alglib.alglibexception(emsg);
1447                sixbits[sixbitsread] = d;
1448                sixbitsread++;
1449                offs++;
1450            }
1451            if( sixbitsread==0 )
1452                throw new alglib.alglibexception(emsg);
1453            for(i=sixbitsread; i<12; i++)
1454                sixbits[i] = 0;
1455            foursixbits2threebytes(sixbits, 0, bytes, 0);
1456            foursixbits2threebytes(sixbits, 4, bytes, 3);
1457            foursixbits2threebytes(sixbits, 8, bytes, 6);
1458            c = (bytes[sizeof(int)-1] & 0x80)!=0 ? (byte)0xFF : (byte)0x00;
1459            for(i=sizeof(int); i<8; i++)
1460                if( bytes[i]!=c )
1461                    throw new alglib.alglibexception(emsg3264);
1462            for(i=0; i<sizeof(int); i++)
1463                _bytes[i] = bytes[i];       
1464            if( !System.BitConverter.IsLittleEndian )
1465                System.Array.Reverse(_bytes);
1466            return System.BitConverter.ToInt32(_bytes,0);
1467        }   
1468       
1469       
1470        /************************************************************************
1471        This function serializes double value into buffer
1472
1473        v           double value to be serialized
1474        buf         buffer, at least 11 characters wide
1475        offs        offset in the buffer
1476       
1477        after return from this function, offs points to the char's past the value
1478        being read.
1479        ************************************************************************/
1480        private static void double2str(double v, char[] buf, ref int offs)
1481        {
1482            int i;
1483            int[] sixbits = new int[12];
1484            byte[] bytes = new byte[9];
1485
1486            //
1487            // handle special quantities
1488            //
1489            if( System.Double.IsNaN(v) )
1490            {
1491                buf[offs+0] = '.';
1492                buf[offs+1] = 'n';
1493                buf[offs+2] = 'a';
1494                buf[offs+3] = 'n';
1495                buf[offs+4] = '_';
1496                buf[offs+5] = '_';
1497                buf[offs+6] = '_';
1498                buf[offs+7] = '_';
1499                buf[offs+8] = '_';
1500                buf[offs+9] = '_';
1501                buf[offs+10] = '_';
1502                offs += SER_ENTRY_LENGTH;
1503                return;
1504            }
1505            if( System.Double.IsPositiveInfinity(v) )
1506            {
1507                buf[offs+0] = '.';
1508                buf[offs+1] = 'p';
1509                buf[offs+2] = 'o';
1510                buf[offs+3] = 's';
1511                buf[offs+4] = 'i';
1512                buf[offs+5] = 'n';
1513                buf[offs+6] = 'f';
1514                buf[offs+7] = '_';
1515                buf[offs+8] = '_';
1516                buf[offs+9] = '_';
1517                buf[offs+10] = '_';
1518                offs += SER_ENTRY_LENGTH;
1519                return;
1520            }
1521            if( System.Double.IsNegativeInfinity(v) )
1522            {
1523                buf[offs+0] = '.';
1524                buf[offs+1] = 'n';
1525                buf[offs+2] = 'e';
1526                buf[offs+3] = 'g';
1527                buf[offs+4] = 'i';
1528                buf[offs+5] = 'n';
1529                buf[offs+6] = 'f';
1530                buf[offs+7] = '_';
1531                buf[offs+8] = '_';
1532                buf[offs+9] = '_';
1533                buf[offs+10] = '_';
1534                offs += SER_ENTRY_LENGTH;
1535                return;
1536            }
1537           
1538            //
1539            // process general case:
1540            // 1. copy v to array of chars
1541            // 2. set 9th byte to zero in order to simplify conversion to six-bit representation
1542            // 3. convert to little endian (if needed)
1543            // 4. convert to six-bit representation
1544            //    (last 12th element of sixbits is always zero, we do not output it)
1545            //
1546            byte[] _bytes = System.BitConverter.GetBytes((double)v);
1547            if( !System.BitConverter.IsLittleEndian )
1548                System.Array.Reverse(_bytes);
1549            for(i=0; i<sizeof(double); i++)
1550                bytes[i] = _bytes[i];
1551            for(i=sizeof(double); i<9; i++)
1552                bytes[i] = 0;
1553            threebytes2foursixbits(bytes, 0, sixbits, 0);
1554            threebytes2foursixbits(bytes, 3, sixbits, 4);
1555            threebytes2foursixbits(bytes, 6, sixbits, 8);
1556            for(i=0; i<SER_ENTRY_LENGTH; i++)
1557                buf[offs+i] = sixbits2char(sixbits[i]);
1558            offs += SER_ENTRY_LENGTH;
1559        }
1560
1561        /************************************************************************
1562        This function unserializes double value from string
1563
1564        buf         buffer which contains value; leading spaces/tabs/newlines are
1565                    ignored, traling spaces/tabs/newlines are treated as  end  of
1566                    the double value.
1567        offs        offset in the buffer
1568       
1569        after return from this function, offs points to the char's past the value
1570        being read.
1571
1572        This function raises an error in case unexpected symbol is found
1573        ************************************************************************/
1574        private static double str2double(char[] buf, ref int offs)
1575        {
1576            string emsg = "ALGLIB: unable to read double value from stream";
1577            int[] sixbits = new int[12];
1578            byte[]  bytes = new byte[9];
1579            byte[] _bytes = new byte[sizeof(double)];
1580            int sixbitsread, i;
1581           
1582           
1583            //
1584            // skip leading spaces
1585            //
1586            while( buf[offs]==' ' || buf[offs]=='\t' || buf[offs]=='\n' || buf[offs]=='\r' )
1587                offs++;
1588           
1589             
1590            //
1591            // Handle special cases
1592            //
1593            if( buf[offs]=='.' )
1594            {
1595                string s = new string(buf, offs, SER_ENTRY_LENGTH);
1596                if( s==".nan_______" )
1597                {
1598                    offs += SER_ENTRY_LENGTH;
1599                    return System.Double.NaN;
1600                }
1601                if( s==".posinf____" )
1602                {
1603                    offs += SER_ENTRY_LENGTH;
1604                    return System.Double.PositiveInfinity;
1605                }
1606                if( s==".neginf____" )
1607                {
1608                    offs += SER_ENTRY_LENGTH;
1609                    return System.Double.NegativeInfinity;
1610                }
1611                throw new alglib.alglibexception(emsg);
1612            }
1613           
1614            //
1615            // General case:
1616            // 1. read and decode six-bit digits
1617            // 2. check that all 11 digits were read
1618            // 3. set last 12th digit to zero (needed for simplicity of conversion)
1619            // 4. convert to 8 bytes
1620            // 5. convert to big endian representation, if needed
1621            //
1622            sixbitsread = 0;
1623            while( buf[offs]!=' ' && buf[offs]!='\t' && buf[offs]!='\n' && buf[offs]!='\r' && buf[offs]!=0 )
1624            {
1625                int d;
1626                d = char2sixbits(buf[offs]);
1627                if( d<0 || sixbitsread>=SER_ENTRY_LENGTH )
1628                    throw new alglib.alglibexception(emsg);
1629                sixbits[sixbitsread] = d;
1630                sixbitsread++;
1631                offs++;
1632            }
1633            if( sixbitsread!=SER_ENTRY_LENGTH )
1634                throw new alglib.alglibexception(emsg);
1635            sixbits[SER_ENTRY_LENGTH] = 0;
1636            foursixbits2threebytes(sixbits, 0, bytes, 0);
1637            foursixbits2threebytes(sixbits, 4, bytes, 3);
1638            foursixbits2threebytes(sixbits, 8, bytes, 6);
1639            for(i=0; i<sizeof(double); i++)
1640                _bytes[i] = bytes[i];       
1641            if( !System.BitConverter.IsLittleEndian )
1642                System.Array.Reverse(_bytes);       
1643            return System.BitConverter.ToDouble(_bytes,0);
1644        }
1645    }
1646   
1647    /*
1648     * Parts of alglib.smp class which are shared with GPL version of ALGLIB
1649     */
1650    public partial class smp
1651    {
1652        #pragma warning disable 420
1653        public const int AE_LOCK_CYCLES = 512;
1654        public const int AE_LOCK_TESTS_BEFORE_YIELD = 16;
1655       
1656        /*
1657         * This variable is used to perform spin-wait loops in a platform-independent manner
1658         * (loops which should work same way on Mono and Microsoft NET). You SHOULD NEVER
1659         * change this field - it must be zero during all program life.
1660         */
1661        public static volatile int never_change_it = 0;
1662       
1663        /*************************************************************************
1664        Lock.
1665
1666        This class provides lightweight spin lock
1667        *************************************************************************/
1668        public class ae_lock
1669        {
1670            public volatile int is_locked;
1671        }
1672
1673        /********************************************************************
1674        Shared pool: data structure used to provide thread-safe access to pool
1675        of temporary variables.
1676        ********************************************************************/
1677        public class sharedpoolentry
1678        {
1679            public apobject obj;
1680            public sharedpoolentry next_entry;
1681        }
1682        public class shared_pool : apobject
1683        {
1684            /* lock object which protects pool */
1685            public ae_lock pool_lock;
1686   
1687            /* seed object (used to create new instances of temporaries) */
1688            public volatile apobject seed_object;
1689           
1690            /*
1691             * list of recycled OBJECTS:
1692             * 1. entries in this list store pointers to recycled objects
1693             * 2. every time we retrieve object, we retrieve first entry from this list,
1694             *    move it to recycled_entries and return its obj field to caller/
1695             */
1696            public volatile sharedpoolentry recycled_objects;
1697           
1698            /*
1699             * list of recycled ENTRIES:
1700             * 1. this list holds entries which are not used to store recycled objects;
1701             *    every time recycled object is retrieved, its entry is moved to this list.
1702             * 2. every time object is recycled, we try to fetch entry for him from this list
1703             *    before allocating it with malloc()
1704             */
1705            public volatile sharedpoolentry recycled_entries;
1706           
1707            /* enumeration pointer, points to current recycled object*/
1708            public volatile sharedpoolentry enumeration_counter;
1709           
1710            /* constructor */
1711            public shared_pool()
1712            {
1713                ae_init_lock(ref pool_lock);
1714            }
1715           
1716            /* initializer - creation of empty pool */
1717            public override void init()
1718            {
1719                seed_object = null;
1720                recycled_objects = null;
1721                recycled_entries = null;
1722                enumeration_counter = null;
1723            }
1724           
1725            /* copy constructor (it is NOT thread-safe) */
1726            public override apobject make_copy()
1727            {
1728                sharedpoolentry ptr, buf;
1729                shared_pool result = new shared_pool();
1730               
1731                /* create lock */
1732                ae_init_lock(ref result.pool_lock);
1733   
1734                /* copy seed object */
1735                if( seed_object!=null )
1736                    result.seed_object = seed_object.make_copy();
1737               
1738                /*
1739                 * copy recycled objects:
1740                 * 1. copy to temporary list (objects are inserted to beginning, order is reversed)
1741                 * 2. copy temporary list to output list (order is restored back to normal)
1742                 */
1743                buf = null;
1744                for(ptr=recycled_objects; ptr!=null; ptr=ptr.next_entry)
1745                {
1746                    sharedpoolentry tmp = new sharedpoolentry();
1747                    tmp.obj =  ptr.obj.make_copy();
1748                    tmp.next_entry = buf;
1749                    buf = tmp;
1750                }
1751                result.recycled_objects = null;
1752                for(ptr=buf; ptr!=null;)
1753                {
1754                    sharedpoolentry next_ptr = ptr.next_entry;
1755                    ptr.next_entry = result.recycled_objects;
1756                    result.recycled_objects = ptr;
1757                    ptr = next_ptr;
1758                }
1759   
1760                /* recycled entries are not copied because they do not store any information */
1761                result.recycled_entries = null;
1762   
1763                /* enumeration counter is reset on copying */
1764                result.enumeration_counter = null;
1765   
1766                return result;
1767            }
1768        }
1769       
1770
1771        /************************************************************************
1772        This function performs given number of spin-wait iterations
1773        ************************************************************************/
1774        public static void ae_spin_wait(int cnt)
1775        {
1776            /*
1777             * these strange operations with ae_never_change_it are necessary to
1778             * prevent compiler optimization of the loop.
1779             */
1780            int i;
1781           
1782            /* very unlikely because no one will wait for such amount of cycles */
1783            if( cnt>0x12345678 )
1784                never_change_it = cnt%10;
1785           
1786            /* spin wait, test condition which will never be true */
1787            for(i=0; i<cnt; i++)
1788                if( never_change_it>0 )
1789                    never_change_it--;
1790        }
1791
1792
1793        /************************************************************************
1794        This function causes the calling thread to relinquish the CPU. The thread
1795        is moved to the end of the queue and some other thread gets to run.
1796        ************************************************************************/
1797        public static void ae_yield()
1798        {
1799            System.Threading.Thread.Sleep(0);
1800        }
1801
1802        /************************************************************************
1803        This function initializes ae_lock structure and sets lock in a free mode.
1804        ************************************************************************/
1805        public static void ae_init_lock(ref ae_lock obj)
1806        {
1807            obj = new ae_lock();
1808            obj.is_locked = 0;
1809        }
1810
1811
1812        /************************************************************************
1813        This function acquires lock. In case lock is busy, we perform several
1814        iterations inside tight loop before trying again.
1815        ************************************************************************/
1816        public static void ae_acquire_lock(ae_lock obj)
1817        {
1818            int cnt = 0;
1819            for(;;)
1820            {
1821                if( System.Threading.Interlocked.CompareExchange(ref obj.is_locked, 1, 0)==0 )
1822                    return;
1823                ae_spin_wait(AE_LOCK_CYCLES);
1824                cnt++;
1825                if( cnt%AE_LOCK_TESTS_BEFORE_YIELD==0 )
1826                    ae_yield();
1827            }
1828        }
1829
1830
1831        /************************************************************************
1832        This function releases lock.
1833        ************************************************************************/
1834        public static void ae_release_lock(ae_lock obj)
1835        {
1836            System.Threading.Interlocked.Exchange(ref obj.is_locked, 0);
1837        }
1838
1839
1840        /************************************************************************
1841        This function frees ae_lock structure.
1842        ************************************************************************/
1843        public static void ae_free_lock(ref ae_lock obj)
1844        {
1845            obj = null;
1846        }
1847       
1848       
1849        /************************************************************************
1850        This function returns True, if internal seed object was set.  It  returns
1851        False for un-seeded pool.
1852
1853        dst                 destination pool (initialized by constructor function)
1854
1855        NOTE: this function is NOT thread-safe. It does not acquire pool lock, so
1856              you should NOT call it when lock can be used by another thread.
1857        ************************************************************************/
1858        public static bool ae_shared_pool_is_initialized(shared_pool dst)
1859        {
1860            return dst.seed_object!=null;
1861        }
1862
1863
1864        /************************************************************************
1865        This function sets internal seed object. All objects owned by the pool
1866        (current seed object, recycled objects) are automatically freed.
1867
1868        dst                 destination pool (initialized by constructor function)
1869        seed_object         new seed object
1870
1871        NOTE: this function is NOT thread-safe. It does not acquire pool lock, so
1872              you should NOT call it when lock can be used by another thread.
1873        ************************************************************************/
1874        public static void ae_shared_pool_set_seed(shared_pool dst, alglib.apobject seed_object)
1875        {
1876            dst.seed_object = seed_object.make_copy();
1877            dst.recycled_objects = null;
1878            dst.enumeration_counter = null;
1879        }
1880
1881
1882        /************************************************************************
1883        This  function  retrieves  a  copy  of  the seed object from the pool and
1884        stores it to target variable.
1885
1886        pool                pool
1887        obj                 target variable
1888       
1889        NOTE: this function IS thread-safe.  It  acquires  pool  lock  during its
1890              operation and can be used simultaneously from several threads.
1891        ************************************************************************/
1892        public static void ae_shared_pool_retrieve<T>(shared_pool pool, ref T obj) where T : alglib.apobject
1893        {
1894            alglib.apobject new_obj;
1895           
1896            /* assert that pool was seeded */
1897            alglib.ap.assert(pool.seed_object!=null, "ALGLIB: shared pool is not seeded, PoolRetrieve() failed");
1898           
1899            /* acquire lock */
1900            ae_acquire_lock(pool.pool_lock);
1901           
1902            /* try to reuse recycled objects */
1903            if( pool.recycled_objects!=null )
1904            {
1905                /* retrieve entry/object from list of recycled objects */
1906                sharedpoolentry result = pool.recycled_objects;
1907                pool.recycled_objects = pool.recycled_objects.next_entry;
1908                new_obj = result.obj;
1909                result.obj = null;
1910               
1911                /* move entry to list of recycled entries */
1912                result.next_entry = pool.recycled_entries;
1913                pool.recycled_entries = result;
1914               
1915                /* release lock */
1916                ae_release_lock(pool.pool_lock);
1917               
1918                /* assign object to smart pointer */
1919                obj = (T)new_obj;
1920               
1921                return;
1922            }
1923               
1924            /*
1925             * release lock; we do not need it anymore because
1926             * copy constructor does not modify source variable.
1927             */
1928            ae_release_lock(pool.pool_lock);
1929           
1930            /* create new object from seed */
1931            new_obj = pool.seed_object.make_copy();
1932               
1933            /* assign object to pointer and return */
1934            obj = (T)new_obj;
1935        }
1936
1937
1938        /************************************************************************
1939        This  function  recycles object owned by the source variable by moving it
1940        to internal storage of the shared pool.
1941
1942        Source  variable  must  own  the  object,  i.e.  be  the only place where
1943        reference  to  object  is  stored.  After  call  to  this function source
1944        variable becomes NULL.
1945
1946        pool                pool
1947        obj                 source variable
1948
1949        NOTE: this function IS thread-safe.  It  acquires  pool  lock  during its
1950              operation and can be used simultaneously from several threads.
1951        ************************************************************************/
1952        public static void ae_shared_pool_recycle<T>(shared_pool pool, ref T obj) where T : alglib.apobject
1953        {
1954            sharedpoolentry new_entry;
1955           
1956            /* assert that pool was seeded */
1957            alglib.ap.assert(pool.seed_object!=null, "ALGLIB: shared pool is not seeded, PoolRecycle() failed");
1958           
1959            /* assert that pointer non-null */
1960            alglib.ap.assert(obj!=null, "ALGLIB: obj in ae_shared_pool_recycle() is NULL");
1961           
1962            /* acquire lock */
1963            ae_acquire_lock(pool.pool_lock);
1964           
1965            /* acquire shared pool entry (reuse one from recycled_entries or malloc new one) */
1966            if( pool.recycled_entries!=null )
1967            {
1968                /* reuse previously allocated entry */
1969                new_entry = pool.recycled_entries;
1970                pool.recycled_entries = new_entry.next_entry;
1971            }
1972            else
1973            {
1974                /*
1975                 * Allocate memory for new entry.
1976                 *
1977                 * NOTE: we release pool lock during allocation because new() may raise
1978                 *       exception and we do not want our pool to be left in the locked state.
1979                 */
1980                ae_release_lock(pool.pool_lock);
1981                new_entry = new sharedpoolentry();
1982                ae_acquire_lock(pool.pool_lock);
1983            }
1984           
1985            /* add object to the list of recycled objects */
1986            new_entry.obj = obj;
1987            new_entry.next_entry = pool.recycled_objects;
1988            pool.recycled_objects = new_entry;
1989           
1990            /* release lock object */
1991            ae_release_lock(pool.pool_lock);
1992           
1993            /* release source pointer */
1994            obj = null;
1995        }
1996
1997
1998        /************************************************************************
1999        This function clears internal list of  recycled  objects,  but  does  not
2000        change seed object managed by the pool.
2001
2002        pool                pool
2003
2004        NOTE: this function is NOT thread-safe. It does not acquire pool lock, so
2005              you should NOT call it when lock can be used by another thread.
2006        ************************************************************************/
2007        public static void ae_shared_pool_clear_recycled(shared_pool pool)
2008        {
2009            pool.recycled_objects = null;
2010        }
2011
2012
2013        /************************************************************************
2014        This function allows to enumerate recycled elements of the  shared  pool.
2015        It stores reference to the first recycled object in the smart pointer.
2016
2017        IMPORTANT:
2018        * in case target variable owns non-NULL value, it is rewritten
2019        * recycled object IS NOT removed from pool
2020        * target variable DOES NOT become owner of the new value; you can use
2021          reference to recycled object, but you do not own it.
2022        * this function IS NOT thread-safe
2023        * you SHOULD NOT modify shared pool during enumeration (although you  can
2024          modify state of the objects retrieved from pool)
2025        * in case there is no recycled objects in the pool, NULL is stored to obj
2026        * in case pool is not seeded, NULL is stored to obj
2027
2028        pool                pool
2029        obj                 reference
2030        ************************************************************************/
2031        public static void ae_shared_pool_first_recycled<T>(shared_pool pool, ref T obj) where T : alglib.apobject
2032        {   
2033            /* modify internal enumeration counter */
2034            pool.enumeration_counter = pool.recycled_objects;
2035           
2036            /* exit on empty list */
2037            if( pool.enumeration_counter==null )
2038            {
2039                obj = null;
2040                return;
2041            }
2042           
2043            /* assign object to smart pointer */
2044            obj = (T)pool.enumeration_counter.obj;
2045        }
2046
2047
2048        /************************************************************************
2049        This function allows to enumerate recycled elements of the  shared  pool.
2050        It stores pointer to the next recycled object in the smart pointer.
2051
2052        IMPORTANT:
2053        * in case target variable owns non-NULL value, it is rewritten
2054        * recycled object IS NOT removed from pool
2055        * target pointer DOES NOT become owner of the new value
2056        * this function IS NOT thread-safe
2057        * you SHOULD NOT modify shared pool during enumeration (although you  can
2058          modify state of the objects retrieved from pool)
2059        * in case there is no recycled objects left in the pool, NULL is stored.
2060        * in case pool is not seeded, NULL is stored.
2061
2062        pool                pool
2063        obj                 target variable
2064        ************************************************************************/
2065        public static void ae_shared_pool_next_recycled<T>(shared_pool pool, ref T obj) where T : alglib.apobject
2066        {   
2067            /* exit on end of list */
2068            if( pool.enumeration_counter==null )
2069            {
2070                obj = null;
2071                return;
2072            }
2073           
2074            /* modify internal enumeration counter */
2075            pool.enumeration_counter = pool.enumeration_counter.next_entry;
2076           
2077            /* exit on empty list */
2078            if( pool.enumeration_counter==null )
2079            {
2080                obj = null;
2081                return;
2082            }
2083           
2084            /* assign object to smart pointer */
2085            obj = (T)pool.enumeration_counter.obj;
2086        }
2087
2088
2089        /************************************************************************
2090        This function clears internal list of recycled objects and  seed  object.
2091        However, pool still can be used (after initialization with another seed).
2092
2093        pool                pool
2094        state               ALGLIB environment state
2095
2096        NOTE: this function is NOT thread-safe. It does not acquire pool lock, so
2097              you should NOT call it when lock can be used by another thread.
2098        ************************************************************************/
2099        public static void ae_shared_pool_reset(shared_pool pool)
2100        {   
2101            pool.seed_object = null;
2102            pool.recycled_objects = null;
2103            pool.enumeration_counter = null;
2104        }
2105    }
2106}
2107public partial class alglib
2108{
2109    public partial class smp
2110    {
2111        public static int cores_count = 1;
2112        public static volatile int cores_to_use = 0;
2113    }
2114    public class smpselftests
2115    {
2116        public static bool runtests()
2117        {
2118            return true;
2119        }
2120    }
2121    public static void setnworkers(int nworkers)
2122    {
2123        alglib.smp.cores_to_use = nworkers;
2124    }
2125}
Note: See TracBrowser for help on using the repository browser.