/SiemensSuite/TotInfo/TotInfo/Version15.cs
C# | 380 lines | 276 code | 71 blank | 33 comment | 40 complexity | d8cce814cecb6152dcf3ca6f87c1346e MD5 | raw file
1using System;
2using System.Collections.Generic;
3using System.Linq;
4using System.Text;
5using System.IO;
6using System.Text.RegularExpressions;
7using Edu.Nlu.Sir.Siemens.Shared;
8namespace Edu.Nlu.Sir.Siemens.TotInfo
9{
10 public class BaseVersion : ITotInfo, FaultyVersion
11 {
12 public int[] FaultLines { get { return new int[] { -1 }; } }
13 public FaultType FaultType { get { return FaultType.NO_FAULT; } }
14 public const int MAXLINE = 256;
15 public const int MAXTBL = 1000;
16 public const char COMMENT = '#';
17 public const int NULL = 0;
18 public const int EXIT_FAILURE = -1;
19 public const int EXIT_SUCCESS = 0;
20
21 public const int ITMAX = 100;
22 public const double EPS = 3.0e-7;
23
24 /// <summary>
25 /// row/column header input line */
26 /// </summary>
27 public char[] line = new char[MAXLINE];
28 /// <summary>
29 /// frequency tallies
30 /// </summary>
31 public long[] f = new long[MAXTBL];
32 /// <summary>
33 /// # of rows
34 /// </summary>
35 public int r;
36 /// <summary>
37 /// # of columns
38 /// </summary>
39 public int c;
40
41
42 void Main(string[] args)
43 {
44 char p; /* input line scan location */
45 int i; /* row index */
46 int j; /* column index */
47 double info; /* computed information measure */
48 int infodf; /* degrees of freedom for information */
49 double totinfo = 0.0; /* accumulated information */
50 int totdf; /* accumulated degrees of freedom */
51
52 totdf = 0;
53 string lineStr;
54 while ((lineStr = Console.ReadLine()) != null) /* start new table */
55 {
56 line = lineStr.ToCharArray();
57 p = line[0];
58 for (int pi = 1; p != '\0' && Char.IsWhiteSpace(p); ++pi)
59 {
60 p = line[pi];
61 }
62
63 if (p == '\0')
64 continue; /* skip blank line */
65
66 if (p == COMMENT)
67 { /* copy comment through */
68 Console.WriteLine(line);
69 continue;
70 }
71
72 string[] tokens = p.ToString().Split(" ".ToCharArray());
73 bool parsable = Int32.TryParse(tokens[0], out r) && Int32.TryParse(tokens[1], out c);
74 if (!parsable)
75 {
76 Console.WriteLine("* invalid row/column line *\n");
77 Environment.Exit(EXIT_FAILURE);
78 }
79
80 if (r * c > MAXTBL)
81 {
82 Console.WriteLine("* table too large *\n");
83 Environment.Exit(EXIT_FAILURE);
84 }
85
86 /* input tallies */
87
88 string valuesString = Console.In.ReadToEnd();
89 string[] valuesLine = valuesString.Split("\n".ToCharArray());
90 string[][] values = (string[][])valuesLine.Select(l => l.Split(" ".ToCharArray()));
91 long value;
92 bool parsed;
93 for (i = 0; i < r; ++i)
94 {
95 for (j = 0; j < c; ++j)
96 {
97 parsed = Int64.TryParse(values[i][j], out value);
98 if (parsed)
99 {
100 f[i * c + j] = value;
101 }
102 else
103 {
104 Console.Write("* EOF in table *\n");
105 Environment.Exit(EXIT_FAILURE);
106 }
107 }
108 }
109
110 /* compute statistic */
111 info = InfoTbl(r, c, f,out infodf);
112
113 /* print results */
114
115 if (info >= 0.0)
116 {
117 Console.Write(String.Format("2info = {0,000.00}\tdf = {1,##}\tq = {2,000.0000}%7.4f\n",
118 info, infodf,
119 QChiSq(info, infodf))
120 );
121 totinfo += info;
122 totdf += infodf;
123 }
124 else
125 Console.Write(info < -3.5 ? "out of memory\n"
126 : info < -2.5 ? "table too small\n"
127 : info < -1.5 ? "negative freq\n"
128 : "table all zeros\n");
129 }
130
131 if (totdf <= 0)
132 {
133 Console.Write("\n*** no information accumulated ***\n");
134 Environment.Exit(EXIT_FAILURE);
135 }
136
137 Console.Write("\ntotal 2info = %5.2f\tdf = %2d\tq = %7.4f\n",
138 totinfo, totdf,
139 QChiSq(totinfo, totdf)
140 );
141 Environment.Exit(EXIT_SUCCESS);
142 }
143
144 public override double LGamma(double x)
145 {
146 double[] cof =
147 {
148 76.18009173, -86.50532033, 24.01409822,
149 -1.231739516, 0.120858003e-2, -0.536382e-5
150 };
151 double tmp, ser;
152 int j;
153
154
155 if (--x < 0.0) /* use reflection formula for accuracy */
156 {
157 double pix = Math.PI * x;
158
159 return Math.Log(pix / Math.Sin(pix)) - LGamma(1.0 - x);
160 }
161
162 tmp = x + 5.5;
163 tmp -= (x + 0.5) * Math.Log(tmp);
164
165 ser = 1.0;
166
167 for (j = 0; j < 6; ++j)
168 ser += cof[j] / ++x;
169
170 return -tmp + Math.Log(2.50662827465 * ser);
171 }
172
173
174
175 public override double
176 gser(double a, double x)
177 {
178 double ap, del, sum;
179 int n;
180
181
182 if (x <= 0.0)
183 return 0.0;
184
185 del = sum = 1.0 / (ap = a);
186
187 for (n = 1; n <= ITMAX; ++n)
188 {
189 sum += del *= x / ++ap;
190
191 if (Math.Abs(del) < Math.Abs(sum) * EPS)
192 return sum * Math.Exp(-x + a * Math.Log(x) - LGamma(a));
193 }
194
195 throw new ApplicationException("Execution should not have reached this line");
196 }
197
198 public override double
199 gcf(double a, double x)
200 {
201 int n;
202 double gold = 0.0, fac = 1.0, b1 = 1.0,
203 b0 = 0.0, a0 = 1.0, a1 = x;
204
205 for (n = 1; n <= ITMAX; ++n)
206 {
207 double anf;
208 double an = (double)n;
209 double ana = an - a;
210
211 a0 = (a1 + a0 * ana) * fac;
212 b0 = (b1 + b0 * ana) * fac;
213 anf = an * fac;
214 b1 = x * b0 + anf * b1;
215 a1 = x * a0 + anf * a1;
216
217 if (a1 != 0.0)
218 { /* renormalize */
219 double g = b1 * (fac = 1.0 / a1);
220
221 gold = g - gold;
222
223 if (Math.Abs(gold) < EPS * Math.Abs(g))
224 return Math.Exp(-x + a * Math.Log(x) - LGamma(a)) * g;
225
226 gold = g;
227 }
228 }
229 throw new ApplicationException("Execution should not have reached this line");
230 }
231
232 public override double
233 QGamma(double a, double x)
234 {
235
236 return x < a + 1.0 ? 1.0 - gser(a, x) : gcf(a, x);
237 }
238
239 public override double
240 QChiSq(double chisq, int df)
241 {
242 return QGamma((double)df / 2.0, chisq / 2.0);
243 }
244
245
246
247
248 /// <summary>
249 /// InfoTbl -- Kullback's information measure for a 2-way contingency table
250 ///last edit: 88/09/19 D A Gwyn
251 ///SCCS ID: @(#)info.c 1.1 (edited for publication)
252 ///Special return values:
253 ///-1.0 entire table consisted of 0 entries
254 ///-2.0 invalid table entry (frequency less than 0)
255 ///-3.0 invalid table dimensions (r or c less than 2)
256 ///-4.0 unable to allocate enough working storage
257 /// </summary>
258 /// <param name="r"># rows in table</param>
259 /// <param name="c"># columns in table </param>
260 /// <param name="f">r*c frequency tallies</param>
261 /// <param name="pdf">return # d.f. for chi-square</param>
262 /// <returns></returns>
263 public override double
264 InfoTbl(int r, int c, long[] f, out int pdf)
265 {
266 int i; /* row index */
267 int j; /* column index */
268 double N; /* (double)n */
269 double info; /* accumulates information measure */
270 double[] xi; /* row sums */
271 double[] xj; /* col sums */
272 int rdf = r - 1; /* row degrees of freedom */
273 int cdf = c - 1; /* column degrees of freedom */
274
275 if (rdf <= 0 || cdf <= 0)
276 {
277 pdf = 0;
278 info = -3.0;
279 goto ret3;
280 }
281
282 pdf = rdf * cdf; /* total degrees of freedom */
283 try {
284 xi = new double[r];
285 }
286 catch(OutOfMemoryException)
287 {
288 info = -4.0;
289 goto ret3;
290 }
291 try {
292 xj = new double[c];
293 }
294 catch(OutOfMemoryException)
295 {
296 info = -4.0;
297 goto ret2;
298 }
299
300 /* compute row sums and total */
301
302 N = 0.0;
303
304 for (i = 0; i < r; ++i)
305 {
306 double sum = 0.0; /* accumulator */
307
308 for (j = 0; j < c; ++j)
309 {
310 long k = f[(i) * c + (j)];
311
312 if (k < 0L)
313 {
314 info = -2.0;
315 goto ret1;
316 }
317
318 sum += (double)k;
319 }
320
321 N += xi[i] = sum;
322 }
323
324 if (N <= 0.0)
325 {
326 info = -1.0;
327 goto ret1;
328 }
329
330 /* compute column sums */
331
332 for (j = 0; j < c; ++j)
333 {
334 double sum = 0.0; /* accumulator */
335
336 for (i = 0; i < r; ++i)
337 sum += (double)f[(i) * c + (j)];
338
339 xj[j] = sum;
340 }
341
342 /* compute information measure (four parts) */
343
344 info = N * Math.Log(N); /* part 1 */
345
346 for (i = 0; i < r; ++i)
347 {
348 double pi = xi[i]; /* row sum */
349
350 if (pi > 0.0)
351 info -= pi * Math.Log(pi); /* part 2 */
352
353 for (j = 0; j < c; ++j)
354 {
355 double pij = (double)f[(i) * c + (j)];
356
357 if (pij > 0.0)
358 info += pij * Math.Log(pij); /* part 3 */
359 }
360 }
361
362 for (j = 0; j < c; ++j)
363 {
364 double pj = xj[j]; /* column sum */
365
366 if (pj > 0.0)
367 info -= pj * Math.Log(pj); /* part 4 */
368 }
369
370 info *= 2.0; /* for comparability with chi-square */
371
372 ret1:
373 xj = null;
374 ret2:
375 xi = null;
376 ret3:
377 return info;
378 }
379 }
380}