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