PageRenderTime 28ms CodeModel.GetById 11ms app.highlight 13ms RepoModel.GetById 1ms app.codeStats 0ms

/SiemensSuite/TotInfo/TotInfo/Version10.cs

https://github.com/lukesandberg/PexFaultLocalization
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}