program genhis(data, histog, genhisp, curve, output); (* generalized histogram program by Gary Stormo modified by: Dr. Thomas D. Schneider National Institutes of Health National Cancer Institute Center for Cancer Research Nanobiology Program Molecular Information Theory Group Frederick, Maryland 21702-1201 toms@ncifcrf.gov permanent email: toms@alum.mit.edu (use only if first address fails) http://www.ccrnp.ncifcrf.gov/~toms/ module libraries required: delman, delmods *) label 1; (* end of program *) const (* begin module version *) version = 1.88; (* of genhis.p 2007 Aug 29 2007 Aug 29, 1.88: qqq 2007 Aug 29, 1.87: curve plot documentation inside file 2007 Aug 29, 1.86: curve plot as a real number, not integer! 2007 Aug 03, 1.85: give 0 anyway when can't compute uncertainty 2007 Jul 31, 1.84: bomb on having only one value in data 2007 Jun 14, 1.83: tolerance limit to exp in gaushist solves bug 2007 Jun 14, 1.82: crash on reading column 6 of rsdata file? error in exponentiation (Result too large) (error #700 at 1af03) 2006 Nov 15, 1.81: provide gaussian data 1.80, 2003 Jul 31: tabs now counted as blanks 1.79, 2002 Nov 1: bug: infinite loop with only one data item! 1.78, 2002 Feb 5: protect against bad non-numerical data 1.77, 2002 Feb 5: upgrade documentation 1.76, 2002 Feb 5: allow blank lines in the file. 1.75, 1999 Dec 3: produce curve file 1.74, 1997 Nov 20: previous changes origin before 1984 jan 6 *) (* end module version *) (* begin module describe.genhis *) (* name genhis: general histogram plotter synopsis genhis(data: in, histog: out, genhisp: in, curve: out, output: out) files data: File of numbers to be histogrammed. Header lines that begin with '*' may be copied from this file to 'histog' or may be skipped. The column from which to read the data may also be specified. See the description of the file 'genhisp' to see how to do this. Once the data region has begun, (that is, there is at least one non '*' line), then lines that begin with '*' are also skipped. Lines that are blank or for which there *is* no data column (because one asked for a column further out than is supplied) are skipped. Data are only separated by spaces or by tabs. Other control characters are not separaters. histog: the output histogram. contains the header lines copied from file 'data', plus data about the numbers (min, max, mean, variance and st. dev.), and the plot. may also contain a standard plot. genhisp: parameter input file. this is used to change any of the parameters from default values. any may be changed and they can be specified in any order. the first character on a line tells what parameter is to be set, the other information sets it. the parameters that can be changed, and their line codes: h - sets header reading; this is followed by two integers, the first specifying the number of lines to copy and the second the number of lines to skip; if the first number is <0 those lines beginning with '*' are copied; default is -1 0. c - sets the column of data that is to be analyzed and plotted; the default is column 1; note: a column is any string of nonblank characters; columns are separated by blanks; p - sets the standard plot; poisson and gaussian plots are available and are specified by following the p by either p or g. NEW: 2006 Nov 15: Capital P or G means provide the curve data in the histog. x - sets the x-axis scale; this is to be followed by either an n or an s, and then a number; if n, then the number of intervals on the x-axis is set; if s, then the size of intervals is set; default is to set the number of intervals to constant 'defslots'. r - sets the range of data to be plotted; this is followed by two numbers which specify the subrange of the data for which the plot is desired; default is to plot all the data. curve: gives the standard plot curve as position and value output: for messages to the user. description This program takes numerical data from a file and plots a histogram of those data. It also calculates the min, max, mean and variance of the data. If desired, the user may get a standard plot, based on the mean and variance, plotted along with the data. The user may specify the size or number of intervals on the x-axis. The y-axis is automatically scaled to fit on a page. The scaling factor is reported to the user. example try file datat7 with genhisp of x n 20 p g see also {Program to plot the graph in postscript:} genpic.p {An x-y plotting program:} xyplo.p author Gary Stormo bugs Try different x axis intervals: regular spikes can be data artifacts! MAJOR BUG: If the data column has identical values, then the program goes into an infinite loop and will generate an infinite sized histog file. Fortunately this is not an interesting histogram ... technical notes The constant 'pageheight' is used to set the scaling factor so that the plots do not exceede the size of a page. As of 1994 Aug 10 the SEM is calculated by dividing the standard deviation by sqrt(n-1). *) (* end module describe.genhis *) (* begin module universal.const *) e = 2.71828459045; pi = 3.14159265358974; (* end module universal.const *) (* begin module histogram.const *) maxslots = 20000; (* maximum number of slots in the histogram *) defslots = 100; (* the default number of slots *) pageheight = 53; (* was 105; *) (* room on page for plotting; used to set scaling. pageheight should be set to 27 less than the number of spaces allowed across a page for your printer. *) dch = '+'; (* character for plotting data *) sch = ':'; (* character for plotting standard *) bch = '*'; (* character for plotting coincidence of standard and data *) wid = 12; (* width of numbers to use for data *) dec = 8; (* number of decimal places to use for data *) (* end module histogram.const *) type (* begin module histogram.type *) plots = (none,gaussian,poisson); (* the types of standard plots that can be done *) endpoints = (start,stop); (* defines the range for plotting *) histarray = record numbers: array[1..maxslots] of integer; (* rounded histogram arrays *) rnumbers: array[1..maxslots] of real; (* histogram as reals *) range: array[endpoints] of real; (* range of the data recorded *) interval: real; (* the size of the histogram slots *) slots: integer; (* the number of slots *) end; (* end module histogram.type *) var (* begin module genhis.var *) genhisp, (* parameter file *) data, (* input file of numbers *) histog, (* output histogram *) curve: text; (* output standard plot *) (* end module genhis.var *) (* begin module histogram.var *) dataarray: histarray; (* the histogram of the data *) stanarray: histarray; (* histogram of standard values *) hlines, (* the number of header lines in the numbers file *) column, (* of the data file to read *) entries: integer; (* the number of numbers in the file *) min, (* minimum of the numbers *) max, (* maximum of the numbers *) ave, (* the average (mean) number *) scaling, (* for plotting on the page *) vari: real; (* the variance of the numbers *) uncertainty: real; (* Shannon's H function for the frequencies *) standplot: plots; (* type of standard plot to use *) splot: boolean; (* true if we plot a standard curve with the data *) curvedata: boolean; (* give curve data *) (* end module histogram.var *) (* begin module package.primitive *) (* ************************************************************************ *) (* begin module halt *) procedure halt; (* stop the program. the procedure performs a goto to the end of the program. you must have a label: label 1; declared, and also the end of the program must have this label: 1: end. examples are in the module libraries. this is the only goto in the delila system. *) begin writeln(output,' program halt.'); goto 1 end; (* end module halt version = 'delmod 6.51 85 apr 17 tds/gds' *) (* begin module copyaline *) procedure copyaline(var fin, fout: text); (* copy a line from file fin to file fout *) begin (* copyaline *) while not eoln(fin) do begin fout^ := fin^; put(fout); get(fin) end; readln(fin); writeln(fout); end; (* copyaline *) (* end module copyaline version = 'delmod 6.51 85 apr 17 tds/gds' *) (* begin module copylines *) function copylines(var fin, fout: text; n: integer): integer; (* copy n lines of file fin to file fout. the actual number of lines copied is returned. *) var index: integer; (* the current line number *) begin (* copylines *) index := 0; while (not eof(fin)) and (index < n) do begin copyaline(fin, fout); index := succ(index) end; copylines := index end; (* copylines *) (* end module copylines version = 'delmod 6.51 85 apr 17 tds/gds' *) (* ************************************************************************ *) (* end module package.primitive version = 'delmod 6.51 85 apr 17 tds/gds' *) (* begin module skipblanks *) (* 2003 July 31: tab is considered a blank character *) function isblank(c: char): boolean; (* is the character c blank or tab? *) begin isblank := (c = ' ') or (c = chr(9)) end; procedure skipblanks(var thefile: text); (* skip over blanks until a non-blank, or end of line, is found *) begin while isblank(thefile^) and not eoln(thefile) do get(thefile); end; procedure skipnonblanks(var thefile: text); (* skip over nonblanks until a blank, or end of line, is found *) begin while (not isblank(thefile^)) and not eoln(thefile) do get(thefile); end; { OLD PROCEDURES: tab previously was not a blank. procedure skipblanks(var thefile: text); (* skip over blanks until a non-blank, or end of line, is found *) begin while (thefile^ = ' ') and not eoln(thefile) do get(thefile); end; procedure skipnonblanks(var thefile: text); (* skip over nonblanks until a blank, or end of line, is found *) begin while (thefile^ <> ' ') and not eoln(thefile) do get(thefile); end; } procedure skipcolumn(var thefile: text); (* skip over a data column *) begin skipblanks(thefile); skipnonblanks(thefile) end; (* end module skipblanks version = 4.59; (@ of prgmod.p 2003 Jul 31 *) (* begin module genhis.setparam *) procedure setparam; (* read from the parameter file the desired parameters. five parameters can be specified, and each specification is designated by the character that begins the line. unspecified parameters are set to default values. parameters: header lines - start line with 'h', followed by two integers; the first integer sets the number of lines to copy from 'numbers' to 'histog' (if <0 then lines beginning with '*' are copied); the second integer specifies the number of lines to skip over in 'numbers'; defaults are -1, 0, respectively. standard plot - start line with 'p'; if this is followed by a 'g' then a gaussian standard will be plotted, otherwise none. 'P' and 'G' mean to provide the data in histog. x-axis scale - start line with 'x' followed by 'n' or 's' followed by an integer; if 'n' then the number of intervals on the x-axis is set by the integer, if 's' then the size of the intervals is set. data column - start line with 'c' followed by an integer which sets the column in the input file that is to be used for the data; the default is 1. range of plot - start line with 'r' followed by two integers; the first integer specifies at what lower limit of data the plotting is to begin, and the second integer specifies at what upper limit it is to end; the default is all the data. *) var ch: char; (* a character for reading parameters *) c, s, (* the number of lines to copy and skip in the header *) i: integer; (* generic integer *) begin reset(data); if eof(data) then begin writeln(output,' data file empty'); halt; end; rewrite(histog); writeln(histog,'* genhis ',version:4:2, ': generalized histogram data from'); (* set the default values for parameters *) hlines := 0; column := 1; standplot := none; dataarray.slots := defslots; (* by setting range[start] > range[stop] we make them get set to the min and max of the data *) dataarray.range[start] := 1; dataarray.range[stop] := 0; (* change parameters according to user requests *) reset(genhisp); while not eof(genhisp) do begin read(genhisp,ch); if not (ch in ['h','p','x','c','r']) then begin writeln(output,'genhisp commands must begin with one of:', ' hpxcr'); halt end else case ch of 'h' : begin readln(genhisp,c,s); if (c < 0) then begin (* we copy lines beginning with '*' to histog *) c := 0; while (data^ = '*') do begin write(histog,'* '); copyaline(data,histog); c := succ(c) end end else begin i := copylines(data,histog,c); if (i < c) then begin writeln(output,' parameter error:', ' header too short'); halt; end; end; for i := 1 to s do readln(data); (* update hlines to the number of lines in the data which have been copied and skipped *) hlines := hlines + c + s; end; 'p' : begin skipblanks(genhisp); readln(genhisp,ch); if not (ch in ['g','p','n','G','P']) then begin writeln(output,'genhisp plot command can accept', ' only g G p P or n commands'); halt end; case ch of 'g': standplot := gaussian; 'G': standplot := gaussian; 'p': standplot := poisson; 'P': standplot := poisson; 'n': standplot := none; end; if (ch in ['G','P']) then curvedata := true else curvedata := false; end; 'x' : with dataarray do begin skipblanks(genhisp); read(genhisp,ch); if (ch = 's') then begin readln(genhisp,interval); if (interval <= 0.0) then begin writeln(output,' error in parameter file: ', ' interval size must be positive'); halt; end; slots := 0 end else if ch = 'n' then begin readln(genhisp,slots); if (slots < 1) then begin writeln(output,' error in parameter file: ', 'number of slots must be positive'); halt end else if (slots > maxslots) then begin writeln(output,' error in parameter file: ', 'number of slots must be less than ',maxslots:1); halt end end else begin writeln(output,' x parameter must be followed', ' by "s" or "n".'); halt end end; 'c' : readln(genhisp,column); 'r' : with dataarray do readln(genhisp,range[start],range[stop]); end; end; (* if we have not copied or skipped any lines in the header, i.e., if hlines = 0, then we copy those lines beginning with '*' *) if (hlines = 0) then begin while (data^ = '*') do begin copyaline(data,histog); hlines := succ(hlines); end; end; writeln(histog,'*'); end; (* end module genhis.setparam *) (* begin module histogram *) (* the following modules contain procedures for calculating and plotting histogram data. *) (* begin module histogram.datahist *) (* end module histogram.datahist *) (* begin module histogram.gaushist *) (* end module histogram.gaushist *) (* begin module histogram.poishist *) (* end module histogram.poishist *) (* begin module histogram.hplot *) (* end module histogram.hplot *) (* end module histogram *) (* begin module histogram.datahist *) procedure datahist(var thefile: text; hlines, column: integer; var entries: integer; var min, max, ave, vari: real; var dataarray: histarray); (* skip hlines in thefile to get to the data. now read the data from the specified 'column' (must be an integer). count the number of entries to the end of the file, determine the min, max, sum and average. now determine the size of the intervals or the number of slots, depending on which has been specified. pass through the data a second time, counting the occurrences for each interval and calculating the variance. this procedure uses procedures skipblanks and skipnonblanks. *) var sum, (* of the entries *) num: real; (* a number read from thefile *) int, (* a particular interval *) i,j: integer; (* indecies *) frequency: real; (* a frequency in a slot *) {BUBBA}bubbaline:integer; (* count of the lines by my friend bubba *) begin bubbaline:=1; (* BUBBA starts counting at line 1 *) reset(thefile); if eof(thefile) then begin writeln(output,' data file is empty'); halt; end; for i := 1 to hlines do begin readln(thefile); bubbaline:=succ(bubbaline); (* BUBBA *) end; entries := 0; sum := 0; (* first pass through the data to calculate min, max and ave *) while not eof(thefile) do begin (* skip lines that begin with '*' *) if thefile^ <> '*' then begin (* get to the proper column for reading the data *) for i := 2 to column do begin skipblanks(thefile); skipnonblanks(thefile); end; skipblanks(thefile); (* 2002 Feb 5: Oh, come on!! Just skip the darn line!!! if eoln(thefile) then begin writeln(output,'parameter error: data column ', column:1,' does not exist'); writeln(output,entries:1,' entries successfully read'); halt; end; *) (* skip blank lines 2002 Feb 5 *) if eoln(thefile) then begin { writeln(output,'skipping blank line'); } readln(thefile); bubbaline:=succ(bubbaline); (* BUBBA *) { writeln(output,'skipped blank line',bubbaline:1); (* BUBBA *) } end else begin { writeln(output,'BUBBA MY OL'' BUDDY!! 2 ',bubbaline:1); (* BUBBA *) writeln(output,'thefile^=',thefile^); } (* check for bad numerical data: *) if not (thefile^ in ['-','+','0','1','2','3','4','5','6','7','8','9']) then begin writeln(output, 'BUBBA found a BAD NUMBER on', ' line ',bubbaline:1, ' column ',column:1); writeln(output, 'The "number" begins with the character ', '"',thefile^,'".'); writeln(output,'BUBBA MY OL'' BUDDY!! Saves the day again!'); halt end; readln(thefile,num); bubbaline:=succ(bubbaline); (* BUBBA *) { writeln(output,'skipped blank line BUBBA',bubbaline:1); (* BUBBA *) } entries := succ(entries); (* min and max are set to the first entry in the data and then are changed as higher and lower values are read *) if (entries = 1) then begin min := num; max := num end else begin if (num < min) then min := num; if (num > max) then max := num; end; sum := sum + num; end; end else begin readln(thefile); bubbaline:=succ(bubbaline); (* BUBBA *) end end; if entries > 0 then begin ave := sum / entries; (* now do the second pass through the data and count the number of data points in each interval *) reset(thefile); for i := 1 to hlines do readln(thefile); vari := 0; with dataarray do begin (* if range[start] > range[stop], that is the signal that we are going to set the range to the whole data *) if (range[start] > range[stop]) then begin range[start] := min; range[stop] := max; end; (* if slots <= 0, that means we have specified the interval size and we determine the number of slots needed from that, else we know the number of slots and set the interval size *) if (slots > 0) then interval := (range[stop] - range[start]) / slots; (* set the exact number of slots needed for the interval size *) if interval = 0 then interval := 1; slots := trunc((range[stop] - range[start]) / interval) + 1; if slots > maxslots then begin writeln(output,' the calculated number of slots (',slots:1, ') will not fit'); writeln(output,' into the maximum number of slots (', 'constant maxslots = ',maxslots:1,')'); halt end; (* clear the numbers *) for i := 1 to slots do numbers[i] := 0; for i := 1 to slots do rnumbers[i] := 0.0; (* read in the data *) { for i := 1 to entries do begin This old method fails because the entries do not correspond directly to the number of lines if * lines are to be skipped! } while not eof(thefile) do begin (* skip lines that begin with '*' *) if thefile^ <> '*' then begin (* get to the proper column for reading the data *) for j := 2 to column do begin skipblanks(thefile); skipnonblanks(thefile); end; (* 2002 Feb 5: skip blank lines *) if eoln(thefile) then begin readln(thefile); end else begin readln(thefile,num); (* if num is within the range we put it into the proper interval *) if (num >= range[start]) and (num <= range[stop]) then begin int := trunc((num - range[start]) / interval) + 1; numbers[int] := succ(numbers[int]); end; vari := vari + (ave - num) * (ave - num); end end else readln(thefile) end; end; vari := vari / (entries - 1); (* calculate the uncertainty of the distribution *) uncertainty := 0.0; for i := 1 to dataarray.slots do begin with dataarray do begin if numbers[i] > 0.0 then begin frequency := numbers[i] / entries; uncertainty := uncertainty - frequency * ln(frequency); end end end; uncertainty := uncertainty / ln(2); end end; (* end module histogram.datahist *) (* begin module histogram.gaushist *) procedure gaushist(entries: integer; ave, vari: real; var stanarray: histarray); (* fill the stanarray with values assuming the data fits a gaussian distribution. the equation used is from bevington, p.r., 'data reduction and error analysis for the physical sciences', p. 44: prob(x) = exp[-1/2((x-ave)/sd)**2]/sd*sqrt(2*pi). where prob(x) is the probability of getting the value x. to determine the expected number of values inside a given interval, we calculate prob(x) for the mid-point of the interval, multiply by the width of the interval and then multiply this total probability by the number of entries in the data. *) const pi = 3.14159265358974; (* used in calculating the curve *) tolerance = -700; (* limit on exponentiation computation. exp(x) where x < -700 can crash: ./genhis: error in exponentiation (Result too large) (error #700 at 1b6bf) *) var sd, (* standard deviation *) d1, (* first denominator, = 2*sd**2 *) d2, (* second denominator, = sd*sqrt(2*pi) *) x, (* the position for which the expectation is calculated *) ex: real; (* probability of getting a value in the interval *) exinner: real; (* ex := exp(exinner) for precomputing exp *) i: integer; (* index *) begin if vari = 0.0 then begin writeln(output,'WARNING: There is no variance to the data'); writeln(output,'so the Gaussian distribution cannot be plotted.'); writeln(output,'No standard plot will be made.'); writeln(histog,'* no standard plot'); splot := false end else begin sd := sqrt(vari); d1 := 2 * vari; d2 := sd * sqrt(2 * pi); { write(output,', sd = ', sd:6:2); write(output,', d1 = ', d1:6:2); write(output,', d2 = ', d2:6:2); write(output,', interval = ', stanarray.interval:6:2); writeln(output); } with stanarray do for i := 0 to (slots - 1) do begin x := range[start] + (i + 0.5) * interval; { write(output,'bubba 1'); write(output,', i = ', i:1); write(output,', x = ', x:6:2); write(output,', ave = ', ave:6:2); writeln(output); write(output, '(x - ave)/d1 = ', (x - ave)/d1 : 6:2); writeln(output); write(output, '(x - ave)/sd = ', (x - ave)/sd : 6:2); writeln(output); write(output, '((-(x - ave) * (x - ave)/d1) / d2 * interval) = ', ((-(x - ave) * (x - ave)/d1) / d2 * interval): 6:2); writeln(output); } exinner := (-(x - ave) * (x - ave)/d1); if exinner < tolerance then ex := 0.0 else ex := exp(exinner) / d2 * interval; { original before 2007 jun 14: ex := exp(-(x - ave) * (x - ave)/d1) / d2 * interval; writeln(output,'bubba 2'); } rnumbers[i+1] := ex * entries; numbers[i+1] := round(rnumbers[i+1]); { original: numbers[i+1] := round(ex * entries); } end; end; end; (* end module histogram.gaushist *) (* begin module histogram.stanplot *) procedure stanplot(var fout: text; stanarray: histarray); (* print the stanarray distribution to file f *) var x:real; (* the position for which the expectation is calculated *) i: integer; (* index *) begin rewrite(fout); writeln(fout,'* position value'); with stanarray do for i := 0 to (slots - 1) do begin x := range[start] + (i + 0.5) * interval; { writeln(fout, x:wid:dec, ' ', numbers[i+1]:wid:dec) gpc says: warning: second field width allowed only when writing `real' type } writeln(fout, x:wid, ' ', rnumbers[i+1]:wid:dec); { writeln(fout, x:wid, ' ', numbers[i+1]:wid); write(fout, x:wid, ' ', numbers[i+1]:wid); if curvedata then write(fout, zzz writeln(fout); } end; end; (* end module histogram.stanplot *) (* begin module histogram.poishist *) procedure poishist(entries: integer; ave: real; var stanarray: histarray); (* fill the stanarray with values assuming the data fits a poisson distribution. the probability distribution is defined as: prob(x) = ave**x * exp(-ave) / 'x-factorial' . note several things: the variance of a poisson distribution is the same as the average; the poisson distribution is only defined for ave > 0 and for x >= 0; the distribution is only really defined for integer x, but we can also calculate it for reals. since for large values of x 'x-factorial' is hard to calculate we use the following approximation based on sterling's approximation for 'x-factorial': ln(prob(x)) = x - ave + x * ln(ave/x) - ln(2*pi*x)/2 . for x = 1 this approximation is only off by the factor 1.08, and allows us to calculate the prob(x) for all real x > 1. for x between 0 and 1 we use the approximation prob(x) = exp(-ave). if ave <= 0 we skip out of the procedure. *) const pi = 3.14159265358974; (* used in calculating the curve *) var x, (* the position for which the expectation is calculated *) lnex, (* the natural log of ex *) ex: real; (* probability of getting a value in the interval around x *) i: integer; (* index *) begin if (ave > 0) then with stanarray do for i := 0 to (slots - 1) do begin (* get the midpoint of the interval *) x := range[start] + (i + 0.5) * interval; (* the poisson is not defined for x < 0, so we set the array to 0 *) if (x < 0) then ex := 0 (* for x between 0 and 1 use ex = exp(-ave) *) else if (x <= 1) then ex := exp(-ave) else begin lnex := x - ave + x * ln(ave/x) - ln(2*pi*x)/2; ex := exp(lnex) * interval; end; { original: numbers[i+1] := round(ex * entries) } rnumbers[i+1] := ex * entries; numbers[i+1] := round(rnumbers[i+1]) end else writeln(output,' warning: poisson not defined for ave <= 0', ', procedure poishist called but not used'); end; (* end module histogram.poishist *) (* begin module histogram.hplot *) procedure hplot(var thefile: text; scaling: real; splot: boolean; dataarray, stanarray: histarray; dch, sch, bch: char); (* plot the data (in dataarray) to 'thefile'. if splot is true, also plot the data in stanarray. the data is plotted using the character 'dch', the standard using 'sch' and 'bch' is used for the coincidence of the two. *) var dval, (* the data value *) sval, (* the standard value *) i, j: integer; (* indicies *) begin write(thefile,'* interval histogram'); if curvedata then write(thefile,' curve'); writeln(thefile); write(thefile,'* beginning value '); if curvedata then write(thefile,' data'); writeln(thefile); writeln(thefile,'*'); for i := 1 to dataarray.slots do begin dval := round(dataarray.numbers[i]*scaling); with dataarray do write(thefile, ' ',(range[start] + (i - 1) * interval):wid:dec, ' ',numbers[i]:wid,' '); { zzz } { write(thefile, ' bubba '); } if curvedata then with stanarray do write(thefile, ' ',rnumbers[i]:wid:dec,' '); {qqq ' ',numbers[i]:wid,' '); } if splot then begin sval := round(stanarray.numbers[i]*scaling); if (dval < sval) then begin for j := 1 to dval do write(thefile,dch); for j := (dval + 1) to (sval - 1) do write(thefile,' '); write(thefile,sch) end else if ((dval = sval) and (dval > 0))then begin for j := 1 to (dval - 1) do write(thefile,dch); write(thefile,bch) end else begin for j := 1 to (sval - 1) do write(thefile,dch); if (sval > 0) then write(thefile,bch); for j := (sval + 1) to dval do write(thefile,dch); end end else for j := 1 to dval do write(thefile,dch); writeln(thefile); end; end; (* end module histogram.hplot *) (* begin module genhis.writeit *) procedure writeit; (* write the information to thefile, and plot the results *) var maxval, (* maximum value in the dataarray *) i: integer; (* an index *) sd: real; (* standard deviation *) begin writeln(histog,'* parameters:'); writeln(histog,'* ',column:wid, ' is the data column used'); writeln(histog,'* ',entries:wid, ' numbers are in the file'); writeln(histog,'* ',min:wid:dec, ' is the minimum number'); writeln(histog,'* ',max:wid:dec, ' is the maximum number'); writeln(histog,'* ',ave:wid:dec, ' is the MEAN'); sd := sqrt(vari); writeln(histog,'* ',sd:wid:dec, ' is the STANDARD DEVIATION'); writeln(histog,'* ',sd/sqrt(entries-1.0):wid:dec, ' is the STANDARD ERROR OF THE MEAN (SEM)'); writeln(histog,'* ',vari:wid:dec, ' is the variance'); writeln(histog,'* ',uncertainty:wid:dec,' is the uncertainty in bits'); if sd > 0 then writeln(histog,'* ',ln(sqrt(2*pi*e)*sd)/ln(2):wid:dec, ' is the computed uncertainty in bits (Shannon p.57)') else writeln(histog,'* 0 [sd <= 0, can''t compute uncertainty from sd]'); writeln(histog,'*'); writeln(histog,'* ',dataarray.range[start]:wid:dec,' to ', dataarray.range[stop]:wid:dec, ' is the range of data plotted'); writeln(histog,'* ',dataarray.interval:wid:dec, ' is the x-axis interval'); writeln(histog,'* ',dataarray.slots:6,' is the number of intervals'); (* find the maximum histogram value and set the scaling *) maxval := 0; for i := 1 to dataarray.slots do if (maxval < dataarray.numbers[i]) then maxval := dataarray.numbers[i]; (* prevent a bomb if no values in the interval *) if maxval = 0 then maxval := pageheight; scaling := pageheight / maxval; if (scaling > 1) then scaling := 1; writeln(histog,'* ',maxval:wid,' is the maximum y-axis value'); writeln(histog,'* ',scaling:wid:dec,' is the y-axis scale'); (* if a standard plot is wanted, generate it now *) writeln(histog,'*'); if (standplot = none) then begin writeln(histog,'* no standard plot'); splot := false end else begin splot := true; (* set the array parameters up the same *) stanarray.range[start] := dataarray.range[start]; stanarray.range[stop] := dataarray.range[stop]; stanarray.slots := dataarray.slots; stanarray.interval := dataarray.interval; if (standplot = gaussian) then begin gaushist(entries,ave,vari,stanarray); stanplot(curve, stanarray); if splot = true then writeln(histog,'* gaussian standard plotted') end else if (standplot = poisson) then begin (* if ave <= 0 we cannot use the poisson *) if (ave > 0) then begin write(histog,'* poisson standard plotted. '); if (stanarray.range[start] < 0) then begin write(histog,'* warning: poisson not defined for x < 0'); writeln(output,' warning: poisson not defined for x < 0'); end; writeln(histog,'*'); poishist(entries,ave,stanarray) end else begin writeln(output,'warning: poisson not defined for ave <= 0,', ' no standard will be plotted'); splot := false; end; end; end; writeln(histog,'*'); if entries > 0 then hplot(histog,scaling,splot,dataarray,stanarray,dch,sch,bch) else writeln(histog, '* (no numbers are in the file,', ' so no histogram is made)'); end; (* end module genhis.writeit *) begin (* genhis *) (* begin module genhis.body *) writeln(output,' genhis ',version:4:2); setparam; datahist(data,hlines,column,entries,min,max,ave,vari,dataarray); writeit; (* end module genhis.body *) 1: end.