Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.6.0/ALGLIB-3.6.0/ap.cs @ 8803

Last change on this file since 8803 was 8803, checked in by mkommend, 12 years ago

#1968: Added [ThreadStatic] to the RNG of ALGLIB and removed lock from random forest algorithm.

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