Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

ap.cs @ 15529

Visit:

Last change on this file since 15529 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	/**************************************************************************
2	ALGLIB 3.11.0 (source code generated 2017-05-11)
3	Copyright (c) Sergey Bochkanov (ALGLIB project).
4
5	>>> SOURCE LICENSE >>>
6	This program is free software; you can redistribute it and/or modify
7	it under the terms of the GNU General Public License as published by
8	the Free Software Foundation (www.fsf.org); either version 2 of the
9	License, or (at your option) any later version.
10
11	This program is distributed in the hope that it will be useful,
12	but WITHOUT ANY WARRANTY; without even the implied warranty of
13	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14	GNU General Public License for more details.
15
16	A copy of the GNU General Public License is available at
17	http://www.fsf.org/licensing/licenses
18	>>> END OF LICENSE >>>
19	**************************************************************************/
20	using System;
21	public 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.xrhs.x-lhs.yrhs.y, lhs.xrhs.y+lhs.yrhs.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 + 1000System.DateTime.Now.Second + 601000*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 wSystem.Math.Sqrt(1+tt);
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.xz.x-z.yz.y, 2z.xz.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 += rows2; / 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	}
2107	public 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.

Context Navigation

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

Download in other formats: