#!/usr/bin/perl -w
#
# Anand Natrajan (www.anandnatrajan.com)
# Tue May 14 16:00:00 EST 2002
#
# Program to analyse results for generating graphs.

use File::Basename;
use File::Path;
use POSIX ":sys_wait_h";
use Getopt::Long;
use Time::Local;
use IO::Handle;
$| = 1;

$MyName = basename($0);
$DefaultFunc = "mean";
$DefaultProj = "oblique";

sub Usage
{
    print STDERR <<EOF;
Usage: $MyName [-help] [-about] [-v[erbose]]
	-xaxis <name> -yaxis <name> [-zaxis <name>]
	[-func <func>] [-proj <proj>]
	[-keep <name>=<value>]+ [-ignore <name>]+
	[-file <file>]+ [-out <file>]+ [-caption <str>]
	[-xtitle <str>] [-ytitle <str>] [-ztitle <str>]
	[-xlog <float>] [-ylog <float>] [-zlog <float>] [-map <file>]+
EOF

	if ($#_ > -1)
	{
		print STDERR <<EOF;

$MyName generates rough-and-ready graphs from large masses of data.

$MyName assumes that data is fed to it in a particular, strict format.
Every data point must be on its own line. Data must be presented as
<name>=<value> pairs, where <name> is some text-only string representing
the name of the field and <value> is some value for that name. Presumably,
values to be plotted are numeric. Any text on the line that does not
conform to the <name>=<value> pattern is ignored. The <name>=<value> pairs
must be separated by whitespace. An example set of data is:
    Results time=53 size=1024 action=read
    Results time=127 size=2048 action=read
    Results time=84 size=1024 action=write
    Results time=156 size=2048 action=write
The valid <name>s in this set are "time", "size" and "action". Data may be
sent in through files or through STDIN. The entire data is analysed according
to suggestions made with command-line options. The output graphs are meant
to be used as initial templates to be tailored for presentations. A graph
for the above data (in a file "file.dat") can be generated like this:
    $MyName -f file.dat -yaxis time -xaxis size
This command says that the input is in "file.dat", the values for the
field name "time" should be on the Y axis and the values for the field
name "size" should be on the X axis. Implicitly, two plots will be created,
one for action=read and the other for action=write. The plots will be output
only as text on the screen.
    $MyName -f file.dat -yaxis time -xaxis size -keep action=read -out out.spc
This command is the same as the previous one except it says to include only
those lines with action=read and output the graph into "out.spc" which is
assumed to be a space-formatted file.

    -help                   : Print this help screen.
    -about                  : Print information about the author.
    -v[erbose]              : Turn verbose mode on.
    -xaxis <name>           : Plot the values for <name> on the X axis.
    -yaxis <name>           : Plot the values for <name> on the Y axis.
    -zaxis <name>           : Plot the values for <name> on the Z axis.
    -func <func>            : Indicate what to do with multiple Y or Z values
                              for a particular X value. Valid values are
                              "mean", "median", "max", "min", "first",
                              "last", "count", "sum".
                              Default: $DefaultFunc.
    -proj <proj>            : Indicate how to project 3D graphs.
                              Valid values are "oblique", "vanishing",
                              "<n1>[:<n2>],<n3>[:<n4>][,<n5>]".
                              Default: $DefaultProj.
    -keep <name>=<value>    : Keep only those lines that have the
                              <name>=<value> pairs. This option may be
                              repeated. If <name> is repeated, the two
                              <value>s are ORed.
    -ignore <name>          : Ignore the field with <name> in the kept lines.
                              This field will not be used to distinguish
                              lines in the input data.
    -file <file>            : File containing input data. May be repeated,
                              in which case, the effect is the same as
                              concatenating all the files. If no file is
                              specified, input is assumed from STDIN.
    -out <file>             : File containing output graph. May be repeated,
                              in which case, output is sent to each. The
                              extension of the file determines its type.
                              Currently-supported extensions are
                              ".txt" for tab-separated text file,
                              ".tbl" for tbl input file,
                              ".spc" for space-formatted file,
                              ".jgr" for jgraph,
                              ".ps" for PostScript (requires jgraph),
                              ".gif" for GIF (requires jgraph, convert),
                              ".jpg" for JPEG (requires jgraph, convert).
                              ".html" for HTML (requires tbl2html).
                              ".xml" for XML.
    -caption <str>          : Use alternative <string> for caption of graph.
    -xtitle <str>           : Use alternative <string> for title of X axis.
    -ytitle <str>           : Use alternative <string> for title of Y axis.
    -ztitle <str>           : Use alternative <string> for title of Z axis.
    -xlog <float>           : Use <float> as base for logarithmic X axis.
    -ylog <float>           : Use <float> as base for logarithmic Y axis.
    -zlog <float>           : Use <float> as base for logarithmic Z axis.
    -map <file>             : Use <file> for mapping labels on X axis. Format:
        X  <value1>  <label1>
        X  <value2>  <label2>
                              May be repeated, in which case, the effect is
                              the same as concatenating all the files.
                              Using a mapping file turns of automatic labels.

Notes:
    $MyName can generate 4-dimensional graphs, where the fourth dimension
    is colour. Changing the colours based on a scale is expected soon. A
    fifth dimension can be plotted over time, but don't hold your breath.
    The -keep option has the side-effect of moving <name>=<value> pairs
    from the legend to the caption. Therefore, use -keep on a pair that you
    know is a useless filter to shorten the legend at the expense of
    lengthening the caption. A better option is to use -ignore on the <name>
    of the field. Of course, you can override the caption with your own text
    with the -caption option.
    Regular expressions can be supplied in the <value> for a pair with
    -keep, but beware! Quote regexps or else the shell may swallow them
    up. Also, '.*' is bad because it'll transcend the whitespace
    boundaries between <name>=<value> pairs. Use '[^\\s]*' instead, where
    "\\s" stands for whitespace such as space or tab, or use -ignore.
    If graphical output is desired, the values for each <name> specified
    for an axis must be numeric.

EOF
	}
    exit(1);
}

sub About
{
    print STDERR <<EOF;
$MyName authored by Anand Natrajan.
URL: www.anandnatrajan.com.

$MyName is freeware. Please use it as you wish. Suggestions, comments and
criticism can be sent to Anand Natrajan.
EOF
	exit(0);
}

# Print only if verbose mode is on.
sub Vprint
{
	foreach my $line (@_) { print STDERR $line if ($opt_v); }
}

# Apply this function to multiple Y or Z values.
sub func
{
	my $listref = shift;
	if ($opt_func eq "sum" or $opt_func eq "mean")
	{
		my $sum = 0;
		foreach my $val (@$listref)
		{
			$sum += $val;
		}
		return $sum if ($opt_func eq "sum");
		return $sum / ($#{ $listref } + 1) if ($opt_func eq "mean");
	}
	elsif ($opt_func eq "max" or $opt_func eq "min")
	{
		my $max = undef;
		my $min = undef;
		foreach my $val (@$listref)
		{
			$max = $val if (!defined($max) or $val > $max);
			$min = $val if (!defined($min) or $val < $min);
		}
		return $max if ($opt_func eq "max");
		return $min if ($opt_func eq "min");
	}
	elsif ($opt_func eq "first" or $opt_func eq "last" or
		$opt_func eq "median" or $opt_func eq "count")
	{
		my $first = $$listref[0];
		my $last = $$listref[-1];
		my $median = $$listref[int(($#{ $listref } + 1) / 2)];
		my $count = $#{ $listref } + 1;
		return $first if ($opt_func eq "first");
		return $last if ($opt_func eq "last");
		return $median if ($opt_func eq "median");
		return $count if ($opt_func eq "count");
	}
	else
	{
		die "Unknown function $opt_func requested.\n";
	}
}

# Use this function to get logarithms in different bases.
sub logb
{
	my $val = shift; my $base = shift;
	return log($val) / log($base);
}

# Use this function to convert degrees to radians.
sub rad
{
	my $pi = 3.1415926535897932;
	my $val = shift;
	return $pi * $val / 180;
}

# Use this function to convert from 3D value to a 2D value. It's not enough
# to just give a 3D value - you also have to give the range for all three
# axes so that the conversion is proportional. Also, note the
# sleight-of-hand that puts the Z values on the Y axis and smears the Y
# values on the X and Z axes.
sub D3toD2
{
	my $xval = shift; my $yval = shift; my $zval = shift;
	my $xmin = shift; my $ymin = shift; my $zmin = shift;
	my $xmax = shift; my $ymax = shift; my $zmax = shift;

	if ($opt_proj =~ /^persp/)
	{
	}
	else
	{
		# Let's assume that projection is specified as a series of numbers
		# like this: n1[:n2],n3[:n4][,n5]. n1 and n3 are delta and theta
		# respectively, i.e., the real angles made by the Y axis with the
		# X and Z axes respectively. If n2 is specified, then along with n1
		# it defines a range for delta, meaning a vanishing line projection,
		# where the line is parallel to the Z axis. If n4 is specified, it's
		# a vanishing line projection, with the line parallel to the X axis.
		# If both n2 and n4 are specified, it's a vanishing point
		# projection. If neither are specified, it's just an oblique
		# projection. If n5 is specified, it defines how much to scale the Y
		# values down.
		my ($deltas, $thetas, $vanish) = split(/,/, $opt_proj);
		my ($delta1, $delta2) = split(/:/, $deltas);
		my ($theta1, $theta2) = split(/:/, $thetas);
		$delta2 = $delta1 if (!defined($delta2));
		$theta2 = $theta1 if (!defined($theta2));
		my $drange = $delta2 - $delta1;
		my $trange = $theta2 - $theta1;
		$vanish = 3 if (!defined($vanish));

		my $ydiff = $yval - $ymin;
		my $yrange = $ymax - $ymin;
		if (defined($opt_ylog))
		{
			# The difference in the two ratios "normalises" the log scale on
			# the Y-axis.
			$ydiff = logb($yval, $opt_ylog) - logb($ymin, $opt_ylog);
			$yrange = logb($ymax, $opt_ylog) - logb($ymin, $opt_ylog);
		}
		# The yscale defines how far along the Y axis the point is.
		my $yscale = $ydiff / $yrange / $vanish;

		my $xdiff = $xval - $xmin;
		my $xrange = $xmax - $xmin;
		if (defined($opt_xlog))
		{
			# The difference in the two ratios "normalises" the log scale on
			# the X-axis.
			$xdiff = logb($xval, $opt_xlog) - logb($xmin, $opt_xlog);
			$xrange = logb($xmax, $opt_xlog) - logb($xmin, $opt_xlog);
		}
		# The dscale indicates which value to choose for delta.
		my $dscale = $xdiff / $xrange * $drange + $delta1;

		my $zdiff = $zval - $zmin;
		my $zrange = $zmax - $zmin;
		if (defined($opt_zlog))
		{
			# The difference in the two ratios "normalises" the log scale on
			# the Z-axis.
			$zdiff = logb($zval, $opt_zlog) - logb($zmin, $opt_zlog);
			$zrange = logb($zmax, $opt_zlog) - logb($zmin, $opt_zlog);
		}
		# The tscale indicates which value to choose for theta.
		my $tscale = $zdiff / $zrange * $trange + $theta1;

		# The xscale is a mapping of the yscale in terms of the X range.
		my $xscale = $yscale * cos(rad($dscale)) * $xrange;
		# The zscale is a mapping of the yscale in terms of the Z range.
		my $zscale = $yscale * cos(rad($tscale)) * $zrange;

		# Add the scale factors to the existing X and Z values.
		my $xnew = $xval + $xscale;
		$xnew = $opt_xlog ** (logb($xval, $opt_xlog) + $xscale)
			if (defined($opt_xlog));   # Scale Y smear to X's log scale.
		my $ynew = $zval + $zscale;
		$ynew = $opt_zlog ** (logb($zval, $opt_zlog) + $zscale)
			if (defined($opt_zlog));   # Scale Y smear to Z's log scale.
		return ($xnew, $ynew);
	}
}

### MAIN ###

&GetOptions("help" => \$opt_help, "about" => \$opt_about,
	"v|verbose" => \$opt_v, "xaxis=s", "yaxis=s", "zaxis=s", "func=s",
	"proj=s", "keep=s@", "ignore=s@", "f|file=s@", "o|out=s@", "caption=s",
	"xtitle=s", "ytitle=s", "ztitle=s", "xlog=f", "ylog=f", "zlog=f",
	"map=s@") or &Usage;
&Usage(1) if ($opt_help);
&About if ($opt_about);
print STDERR "You must specify -xaxis and -yaxis options.\n\n" and &Usage
	if (!defined($opt_xaxis) or !defined($opt_yaxis));
print STDERR "Base for a log scale must be greater than 0, but not 1.\n\n"
	and &Usage
	if ((defined($opt_xlog) and ($opt_xlog <= 0 or $opt_xlog == 1)) or
		(defined($opt_ylog) and ($opt_ylog <= 0 or $opt_ylog == 1)) or
		(defined($opt_zlog) and ($opt_zlog <= 0 or $opt_zlog == 1)));

# First get in all the input, whether from file or STDIN.
@AllData = ();
if (defined(@opt_f))
{
	foreach my $file (@opt_f)
	{
		open(IN, $file) or die "Couldn't open $file: $!\n";
		push(@AllData, <IN>);
		close IN;
		&Vprint("\nDone reading in file $file.\n");
	}
}
else
{
	push (@AllData, <STDIN>);
	&Vprint("\nDone reading in STDIN.\n");
}
&Vprint("Input data has " . ($#AllData + 1) . " lines.\n");

# Next, parse and swallow the map files, if specified.
%Maps = ();
{ my %dummy = (); $Maps{"X"} = { %dummy }; }
{ my %dummy = (); $Maps{"Y"} = { %dummy }; }
{ my %dummy = (); $Maps{"Z"} = { %dummy }; }
if (defined(@opt_map))
{
	foreach my $mapfile (@opt_map)
	{
		open(MAP, $mapfile) or die "Couldn't open $mapfile: $!\n";
		my @AllMaps = <MAP>;
		close MAP;
		&Vprint("\nDone reading in map $mapfile.\n");
		&Vprint("Map $mapfile has " . ($#AllMaps + 1) . " lines.\n");
		# Okay, now parse the maps and figure out what they mean.
		chomp(@AllMaps);
		foreach my $line (@AllMaps)
		{
			$line =~ s/^\s+//; $line =~ s/\s+$//;
			$line =~ m/([^\s]+)\s+([^\s]+)\s+(.*)/;
			$keyword = $1;
			$value = $2;
			$label = $3;
			if ($keyword !~ /^[xXyYzZ]$/)
			{
				warn "Keyword $keyword unrecognised in $mapfile.\n";
				next;
			}
			$keyword =~ tr/xyz/XYZ/;  # Get them all in uppercase.
			$Maps{$keyword}{$value} = $label;
		}
	}
}

# Back to data. Filter out the lines that are requested only. Add the
# filtering field names to the caption for the graph.
$Caption = "";
@FilterData = @AllData;
&Vprint("Filtered data has " . ($#FilterData + 1) . " lines.\n");
if (defined(@opt_keep))
{
	# Perform processing to ensure that if a name is repeated, its values
	# are ORed.
	my %names = ();
	foreach my $keep (@opt_keep)
	{
		my ($name, $value) = split(/=/, $keep);
		if (defined($names{$name}))
		{
			$names{$name} =~ s/[()]*//g;
			$names{$name} = "($names{$name}|$value)";
		}
		else
		{
			$names{$name} = $value;
		}
	}
	foreach my $name (keys %names)
	{
		&Vprint("\tFiltering for values $names{$name} for name $name.\n");
		@FilterData = grep(/\b$name=$names{$name}\b/, @FilterData);
		$Caption .= " $name=$names{$name}";
		&Vprint("Filtered data has " . ($#FilterData + 1) . " lines.\n");
	}
	# Take out all the pairs that have been used for filtering - their
	# values should be the same, so they can be combined away.
	foreach my $name (keys %names)
	{
		&Vprint("\tFiltering for values $names{$name} for name $name.\n");
		foreach my $line (@FilterData)
		{
			$line =~ s/\b$name=$names{$name}\b//;
		}
	}
}
chomp(@FilterData);
# Remove leading and trailing white spaces in the caption.
$Caption =~ s/^\s+//;
$Caption =~ s/\s+$//;

# After the filtering, from each line remove the fields the user asked to
# ignore. These fields are presumably not data distinguishers.
if (defined(@opt_ignore))
{
	foreach my $ignore (@opt_ignore)
	{
		&Vprint("\tIgnoring for name $ignore.\n");
		foreach my $line (@FilterData)
		{
			# Cull out each field which we're supposed to ignore.
			$line =~ s/\s*\b$ignore=[^\s]*\b//;
		}
	}
}

# After the ignoring, isolate axes values. The remaining <name>s, are
# assumed to refer to separate plots. And the substrings without an '=' are
# discarded.
%Plots = ();
foreach my $line (@FilterData)
{
	my @key = ();
	my $xval = undef;
	my $yval = undef;
	my $zval = undef;
	foreach my $field (split(/\s+/, $line))
	{
		next if ($field !~ /=/);
		my ($name, $value) = split(/=/, $field);
		if ($name eq $opt_xaxis)
		{
			$xval = $value;
		}
		elsif ($name eq $opt_yaxis)
		{
			$yval = $value;
		}
		elsif (defined($opt_zaxis) and $name eq $opt_zaxis)  # 3D
		{
			$zval = $value;
		}
		else
		{
			push(@key, "$name=$value");
		}
	}
	die "Value for X axis variable $opt_xaxis not found in line\n$line\n"
		if (!defined($xval));
	die "Value for Y axis variable $opt_yaxis not found in line\n$line\n"
		if (!defined($yval));
	die "Value for Z axis variable $opt_zaxis not found in line\n$line\n"
		if (defined($opt_zaxis) and !defined($zval));
	my $finalkey = join(" ", sort @key);
	if (!defined($Plots{$finalkey}))
	{
		my @array = ();
		$Plots{$finalkey} = [ @array ];
	}
	my @point = ($xval, $yval, $zval);
	push (@{ $Plots{$finalkey} }, \@point);
}

# Apply the chosen function on multiple Y or Z values. Store the result
# back in the same hash.
$opt_func = $DefaultFunc if (!defined($opt_func));
$opt_proj = $DefaultProj if (!defined($opt_proj));
foreach my $plot (sort keys %Plots)
{
	my %pointhash = ();
	&Vprint("\nKey is $plot\n");
	foreach my $pointref (@{ $Plots{$plot} })
	{
		my $xval = $pointref->[0];
		my $yval = $pointref->[1];
		my $zval = $pointref->[2] if (defined($opt_zaxis));  # 3D
		&Vprint("\tOriginal point is [$xval $yval"
			. (defined($opt_zaxis) ? " $zval" : "") . "]\n");
		# For each point, construct a complex hash to see if the
		# independent variables happen to have multiple values. If they
		# do, then we apply some function on those multiple values.
		if (!defined($opt_zaxis))
		{
			if (!defined($pointhash{$xval}))
			{
				my @ylist = ();
				$pointhash{$xval} = [ @ylist ];
			}
			push(@{ $pointhash{$xval} }, $yval);
		}
		else  # 3D
		{
			if (!defined($pointhash{$xval}))
			{
				my %yhash = ();
				$pointhash{$xval} = { %yhash };
			}
			if (!defined($pointhash{$xval}{$yval}))
			{
				my @zlist = ();
				$pointhash{$xval}{$yval} = [ @zlist ];
			}
			push(@{ $pointhash{$xval}{$yval} }, $zval);
		}
	}
	&Vprint("\tApplying function $opt_func.\n");
	foreach my $xval (sort keys %pointhash)
	{
		if (!defined($opt_zaxis))
		{
			my $result = &func($pointhash{$xval});
			delete $pointhash{$xval};
			$pointhash{$xval} = $result;
			&Vprint("\tNew point is [$xval $result]\n");
		}
		else  # 3D
		{
			foreach my $yval (sort keys %{ $pointhash{$xval} })
			{
				my $result = &func($pointhash{$xval}{$yval});
				delete $pointhash{$xval}{$yval};
				$pointhash{$xval}{$yval} = $result;
				&Vprint("\tNew point is [$xval $yval $result]\n");
			}
		}
	}
	delete $Plots{$plot};
	$Plots{$plot} = { %pointhash };
}

# Now, print the plots.
print "\nGraph Results\n";
print "X Axis: " . (defined($opt_xtitle)
	? "$opt_xtitle ($opt_xaxis)\n" : "$opt_xaxis\n");
print "Y Axis: " . (defined($opt_ytitle)
	? "$opt_ytitle ($opt_yaxis)\n" : "$opt_yaxis\n");
if (!defined($opt_zaxis))
{
	print "Y values function: $opt_func\n";
}
else  # 3D
{
	print "Z Axis: " . (defined($opt_ztitle)
		? "$opt_ztitle ($opt_zaxis)\n" : "$opt_zaxis\n");
	print "Z values function: $opt_func\n";
	print "Z values projection: $opt_proj\n";
}
if ($opt_proj =~ /^obliq/)
{
	$opt_proj = "45:45,45:45";
}
elsif ($opt_proj =~ /^vanis/)
{
	$opt_proj = "45:135,45:135";
}
print "Caption: " . (defined($opt_caption)
	? "$opt_caption ($Caption)\n" : "$Caption\n");
foreach my $plot (sort keys %Plots)
{
	print "\tLegend: $plot\n";
	foreach my $xval (sort {$a <=> $b} keys %{ $Plots{$plot} })
	{
		if (!defined($opt_zaxis))
		{
			print "\t\t$xval, $Plots{$plot}{$xval}\n";
		}
		else  # 3D
		{
			foreach my $yval (sort {$a <=> $b} keys %{ $Plots{$plot}{$xval} })
			{
				print "\t\t$xval, $yval, $Plots{$plot}{$xval}{$yval}\n";
			}
		}
	}
}

# Create the required output files.
&Vprint("Writing " . ($#opt_o + 1) . " output files.\n");
if (defined(@opt_o))
{
	foreach my $outfile (@opt_o)
	{
		&Vprint("Writing output file $outfile.\n");
		# Tab-separated txt file, tbl input file, space-formatted file.
		if ($outfile =~ /\.txt$/ or $outfile =~ /\.tbl$/ or
			$outfile =~ /\.spc$/ or $outfile =~ /\.s?html?/)
		{
			# Save the old file name so we can post-process if necessary.
			my $realfile = $outfile;
			$outfile =~ s/\.spc$/.tbl/;
			$outfile =~ s/\.s?htm?l$/.tbl/;

			my $finalxtitle = "" . (defined($opt_xtitle)
				? $opt_xtitle : $opt_xaxis);
			my $finalytitle = "" . (defined($opt_ytitle)
				? $opt_ytitle : $opt_yaxis);
			my $finalztitle = "" . (defined($opt_ztitle)
				? $opt_ztitle : $opt_zaxis)
				if (defined($opt_zaxis));
			my $finalcaption = "" . (defined($opt_caption)
				? $opt_caption : $Caption);

			open(OUT, ">$outfile") or
				die "Couldn't open $outfile for graph: $!\n";
			# Write the preamble for tbl files only.
			if ($outfile =~ /\.tbl$/)
			{
				if (!defined($opt_zaxis))
				{
					print OUT <<EOF;
.TS
c s s
c c c
l n n.
EOF
				}
				else  # 3D
				{
					print OUT <<EOF;
.TS
c s s s
c c c c
l n n n.
EOF
				}
			}
			print OUT "$finalcaption\nLegend\t$finalxtitle\t$finalytitle"
				. (defined($opt_zaxis) ? "\t$finalztitle" : "") . "\n";

			foreach my $plot (sort keys %Plots)
			{
				foreach my $xval (sort {$a <=> $b} keys %{ $Plots{$plot} })
				{
					if (!defined($opt_zaxis))
					{
						print OUT "$plot\t$xval\t$Plots{$plot}{$xval}\n";
					}
					else  # 3D
					{
						foreach my $yval (sort {$a <=> $b} keys %{ $Plots{$plot}{$xval} })
						{
							print OUT "$plot\t$xval\t$yval\t$Plots{$plot}{$xval}{$yval}\n";
						}
					}
				}
			}
			print OUT ".TE\n"
				if ($outfile =~ /\.tbl$/);
			close OUT;

			# Post-process if necessary to get the proper file.
			if ($realfile =~ /\.spc$/)
			{
				system("tbl < $outfile | nroff -ms | uniq > $realfile");
			}
			elsif ($realfile =~ /\.s?html?$/)
			{
				system("tbl2html < $outfile > $realfile 2> /dev/null");
			}
		}
		# XML file.
		elsif ($outfile =~ /\.xml$/)
		{
			# Save the old file name so we can post-process if necessary.
			my $realfile = $outfile;

			my $finalxtitle = "" . (defined($opt_xtitle)
				? $opt_xtitle : $opt_xaxis);
			my $finalytitle = "" . (defined($opt_ytitle)
				? $opt_ytitle : $opt_yaxis);
			my $finalztitle = "" . (defined($opt_ztitle)
				? $opt_ztitle : $opt_zaxis)
				if (defined($opt_zaxis));
			my $finalcaption = "" . (defined($opt_caption)
				? $opt_caption : $Caption);

			open(OUT, ">$outfile") or
				die "Couldn't open $outfile for graph: $!\n";

			print OUT <<EOF;
<?xml version="1.0" encoding="UTF-8"?>
<graph>
	<title>$finalcaption</title>
	<xaxis>
		<title>$finalxtitle</title>
EOF
			print OUT <<EOF
		<log base="$opt_xlog"/>
EOF
				if (defined($opt_xlog));
			print OUT <<EOF;
	</xaxis>
	<yaxis>
EOF
			print OUT <<EOF
		<log base="$opt_ylog"/>
EOF
				if (defined($opt_ylog));
			print OUT <<EOF;
		<title>$finalytitle</title>
	</yaxis>
EOF

			if (defined($opt_zaxis))
			{
				print OUT <<EOF;
	<zaxis>
		<title>$finalztitle</title>
EOF
				print OUT <<EOF
		<log base="$opt_zlog"/>
EOF
					if (defined($opt_zlog));
				print OUT <<EOF;
	</zaxis>
EOF
			}

			foreach my $plot (sort keys %Plots)
			{
				print OUT <<EOF;
	<plot>
		<title>$plot</title>
EOF
				foreach my $xval (sort {$a <=> $b} keys %{ $Plots{$plot} })
				{
					if (!defined($opt_zaxis))
					{
						print OUT <<EOF;
		<point>
			<x>$xval</x>
			<y>$Plots{$plot}{$xval}</y>
		</point>
EOF
					}
					else  # 3D
					{
						foreach my $yval (sort {$a <=> $b} keys %{ $Plots{$plot}{$xval} })
						{
							print OUT <<EOF;
		<point>
			<x>$xval</x>
			<y>$yval</y>
			<z>$Plots{$plot}{$xval}{$yval}</z>
		</point>
EOF
						}
					}
				}
				print OUT <<EOF;
	</plot>
EOF
			}
			print OUT <<EOF;
</graph>
EOF
			close OUT;
		}
		# jgraph file, ps file, gif file, jpg file.
		elsif ($outfile =~ /\.jgr$/ or $outfile =~ /\.ps$/ or
			$outfile =~ /\.gif$/ or $outfile =~ /\.jpg$/)
		{
			# Save the old file name so we can post-process if necessary.
			my $realfile = $outfile;
			$outfile =~ s/\.ps$/.jgr/;
			$outfile =~ s/\.gif$/.jgr/;
			$outfile =~ s/\.jpg$/.jgr/;

			# Set the thickness of lines for all plots.
			my $linethickness = 1;

			# Set the different line types.
			my @linetypes =
				("solid", "dotted", "dashed", "longdash", "dotdash");
			# my $linetypeindex = rand($#linetypes + 1);
			my $linetypeindex = 0;

			# Set the different markers for points.
			my @marktypes =
				("x", "circle", "box", "cross", "diamond", "triangle", "ellipse");
			# my $marktypeindex = rand($#marktypes + 1);
			my $marktypeindex = 0;

			# Set the different colours for lines as RGB triplets. The
			# choice of colours reflects what is visible and looks good.
			my @colours =
				("0 0 1", "1 0 0", "0 1 0", "1 0 1", "1 1 0", "0 0 0");
			# my $colourindex = rand($#colours + 1);
			my $colourindex = 0;

			# Note that the numbers of linetypes, marktypes and colours
			# are different. This way, as we cycle through them, we can get
			# 5 * 7 * 6 = 210 combinations, which is plenty for one graph.

			my $finalxtitle = "" . (defined($opt_xtitle)
				? "$opt_xtitle\n(* $opt_xaxis *)" : $opt_xaxis);
			my $finalytitle = "" . (defined($opt_ytitle)
				? "$opt_ytitle\n(* $opt_yaxis *)" : $opt_yaxis);
			my $finalztitle = "" . (defined($opt_ztitle)
				? "$opt_ztitle\n(* $opt_zaxis *)" : $opt_zaxis)
				if (defined($opt_zaxis));
			my $finalcaption = "" . (defined($opt_caption)
				? "$opt_caption\n(* $Caption *)" : $Caption);

			# The first thing to do is find the minimum and maximum values
			# along each axis. These values are necessary for setting the
			# limits on the axes as well as scaling, in case we show a 3D
			# graph.
			my $xmin = undef; my $xmax = undef;
			my $ymin = undef; my $ymax = undef;
			my $zmin = undef; my $zmax = undef;
			foreach my $plot (keys %Plots)
			{
				foreach my $xval (keys %{ $Plots{$plot} })
				{
					$xmin = $xval if (!defined($xmin) or $xval < $xmin);
					$xmax = $xval if (!defined($xmax) or $xval > $xmax);
					if (!defined($opt_zaxis))
					{
						my $yval = $Plots{$plot}{$xval};
						$ymin = $yval if (!defined($ymin) or $yval < $ymin);
						$ymax = $yval if (!defined($ymax) or $yval > $ymax);
					}
					else  # 3D
					{
						foreach my $yval (keys %{ $Plots{$plot}{$xval} })
						{
							$ymin = $yval if (!defined($ymin) or $yval < $ymin);
							$ymax = $yval if (!defined($ymax) or $yval > $ymax);
							my $zval = $Plots{$plot}{$xval}{$yval};
							$zmin = $zval if (!defined($zmin) or $zval < $zmin);
							$zmax = $zval if (!defined($zmax) or $zval > $zmax);
						}
					}
				}
			}
			&Vprint("Minimum X value is $xmin.\n");
			&Vprint("Maximum X value is $xmax.\n");
			die "Minimum value cannot be zero or less for log scale.\n"
				if (defined($opt_xlog) and $xmin <= 0);
			&Vprint("Minimum Y value is $ymin.\n");
			&Vprint("Maximum Y value is $ymax.\n");
			die "Minimum value cannot be zero or less for log scale.\n"
				if (defined($opt_ylog) and $ymin <= 0);
			my $xdraw_at = $ymin;
			if (defined($opt_zaxis))
			{
				&Vprint("Minimum Z value is $zmin.\n");
				&Vprint("Maximum Z value is $zmax.\n");
				die "Minimum value cannot be zero or less for log scale.\n"
					if (defined($opt_zlog) and $zmin <= 0);
				$xdraw_at = $zmin;
			}

			# Note the draw_at parameter uses the Z minimum because we plot
			# the Z values on the Y axis.
			my $logstr = "";
			$logstr = "log log_base $opt_xlog" if (defined($opt_xlog));
			open(OUT, ">$outfile") or
				die "Couldn't open $outfile for graph: $!\n";
			print OUT <<EOF;
newgraph
x_translate 0 y_translate 0
title : $finalcaption

(* X axis, with X values *)
xaxis size 5 $logstr draw_at $xdraw_at label : $finalxtitle
(* min $xmin max $xmax *)
no_auto_hash_labels
hash_labels vjc hjr rotate 90
EOF

			# Get the hash labels exactly once for each X value.
			my %hashlabels = ();
			my $foundlabels = 0;
			if (defined(@opt_map))
			{
				my @labels = keys %{ $Maps{"X"} };
				if ($#labels >= 0)
				{
					%hashlabels = %{ $Maps{"X"} };
					$foundlabels = 1;
				}
			}
			if (!$foundlabels)
			{
				foreach my $plot (keys %Plots)
				{
					foreach my $xval (keys %{ $Plots{$plot} })
					{
						$hashlabels{$xval} = $xval;
					}
				}
			}
			foreach my $hash (sort keys %hashlabels)
			{
				print OUT <<EOF;
hash_label at $hash : $hashlabels{$hash}
EOF
			}

			if (!defined($opt_zaxis))
			{
				my $logstr = "";
				$logstr = "log log_base $opt_ylog" if (defined($opt_ylog));
				print OUT <<EOF;

(* Y axis, with Y values *)
yaxis size 5 $logstr draw_at $xmin label : $finalytitle
(* min $ymin max $ymax *)
EOF
				foreach my $plot (sort keys %Plots)
				{
					# Set the line, mark and colour for this curve.
					my $linetype = $linetypes[$linetypeindex];
					$linetypeindex = ($linetypeindex + 1) % ($#linetypes + 1);
					my $marktype = $marktypes[$marktypeindex];
					$marktypeindex = ($marktypeindex + 1) % ($#marktypes + 1);
					my $colour = $colours[$colourindex];
					$colourindex = ($colourindex + 1) % ($#colours + 1);

					print OUT <<EOF;

(* Curve for plot $plot *)
newcurve linethickness $linethickness linetype $linetype
	marktype $marktype color $colour label : $plot
	pts
EOF
					foreach my $xval (sort {$a <=> $b} keys %{ $Plots{$plot} })
					{
						print OUT "\t\t$xval $Plots{$plot}{$xval}\n";
					}
				}
			}
			else  # 3D
			{
				# Once again, note the use of the Z values for the Y axis.
				my $logstr = "";
				$logstr = "log log_base $opt_zlog" if (defined($opt_zlog));
				print OUT <<EOF;

(* Y axis, with Z values *)
yaxis size 5 $logstr draw_at $xmin label : $finalztitle
(* min $zmin max $zmax *)
no_auto_hash_labels
EOF
				# Get the hash labels exactly once for each value. We'll
				# use these for ticks on the artificial Y axis.
				my %xvalues = (); my %yvalues = (); my %zvalues = ();
				foreach my $plot (keys %Plots)
				{
					foreach my $xval (keys %{ $Plots{$plot} })
					{
						$xvalues{$xval} = $xval;
						foreach my $yval (keys %{ $Plots{$plot}{$xval} })
						{
							$yvalues{$yval} = $yval;
							my $zval = $Plots{$plot}{$xval}{$yval};
							$zvalues{$zval} = $zval;
						}
					}
				}

				# Prune the Z values or else the hashes look messy.
				my $sentinel = $zmin;
				my $tithe = ($zmax - $zmin) / 10;
				$tithe = (logb($zmax, $opt_zlog) - logb($zmin, $opt_zlog)) / 10
					if (defined($opt_zlog));
				foreach my $zval (sort {$a <=> $b} keys %zvalues)
				{
					if (!defined($opt_zlog))
					{
						if ($zval >= $sentinel + $tithe and
							$zval <= $zmax - $tithe)
						{
							$sentinel = $zval;
						}
						else
						{
							delete $zvalues{$zval};
						}
					}
					else
					{
						if (logb($zval, $opt_zlog) >= logb($sentinel, $opt_zlog) + $tithe and
							logb($zval, $opt_zlog) <= logb($zmax, $opt_zlog) - $tithe)
						{
							$sentinel = $zval;
						}
						else
						{
							delete $zvalues{$zval};
						}
					}
				}
				$zvalues{$zmin} = $zmin;
				$zvalues{$zmax} = $zmax;

				# Also, for the Z values, which are plotted on the Y axis,
				# put the hash labels where we want them only.
				foreach my $zval (sort keys %zvalues)
				{
					print OUT <<EOF;
hash_label at $zval : $zval
EOF
				}

				# Here's the artificial Y axis, with the hashlabels marked.
				my $ycurve = "";
				print OUT <<EOF;

(* Intermediate axis, with Y values *)
newcurve linethickness 1 linetype solid marktype cross color 0 0 0
	pts
EOF
				my $dist = 0.05;  # Empirical choice for distance of labels.
				# The variable below calculates where to put the labels by
				# simply displacing the X value for the marker.
				foreach my $yval (sort keys %yvalues)
				{
					# Note that we're using the min values of X and Z to
					# get the Y axis.
					my ($xreal, $yreal)
						= &D3toD2($xmax, $yval, $zmin,
							$xmin, $ymin, $zmin, $xmax, $ymax, $zmax);
					my $deltax = $xreal * $dist;
					$deltax = $dist * 4
						if (defined($opt_xlog));
					print OUT "\t\t$xreal $yreal\n";
					$ycurve .= <<EOF;

(* Labels for intermediate axis *)
newcurve marktype text x $deltax y 0 : $yvalues{$yval}
	pts $xreal $yreal
EOF
				}
				# Nice touch, with the Y hashes labelled neatly.
				print OUT $ycurve;
				# Also, let's label the artificial Y axis. The mid point is
				# the linear or logarithmic average.
				my $ymid = $ymin + ($ymax - $ymin) / 2;
				$ymid = $opt_ylog ** (logb($ymin, $opt_ylog)
					+ (logb($ymax, $opt_ylog) - logb($ymin, $opt_ylog)) / 2)
					if (defined($opt_ylog));
				my ($xreal, $yreal)
					= &D3toD2($xmax, $ymid, $zmin,
						$xmin, $ymin, $zmin, $xmax, $ymax, $zmax);
				$dist = .1;  # Empirical choice for distance of title.
				# The variable below once again shows the place to put the
				# title of the fake Y axis. A bit complex for log scale.
				my $deltax = $xreal * $dist;
				$deltax = $dist * 4
					if (defined($opt_xlog));
				print OUT <<EOF;

(* Title for intermediate axis *)
newcurve marktype text vjb hjl rotate 0 fontsize 10 x $deltax y 0 : $finalytitle
	pts $xreal $yreal
EOF

				# Another nice touch is to put a light grey grid on the XY
				# plane so that it's easy to track the points.
				foreach my $xval (sort {$a <=> $b} keys %xvalues)
				{
					next if ($xval == $xmin);
					my ($x1, $y1)
						= &D3toD2($xval, $ymin, $zmin,
							$xmin, $ymin, $zmin, $xmax, $ymax, $zmax);
					my ($x2, $y2)
						= &D3toD2($xval, $ymax, $zmin,
							$xmin, $ymin, $zmin, $xmax, $ymax, $zmax);
					print OUT <<EOF;

(* Grey XY grid line, vertical *)
newcurve linethickness .1 linetype solid marktype none gray 0.1
	pts $x1 $y1 $x2 $y2
EOF
				}
				foreach my $yval (sort {$a <=> $b} keys %yvalues)
				{
					next if ($yval == $ymin);
					my ($x1, $y1)
						= &D3toD2($xmin, $yval, $zmin,
							$xmin, $ymin, $zmin, $xmax, $ymax, $zmax);
					my ($x2, $y2)
						= &D3toD2($xmax, $yval, $zmin,
							$xmin, $ymin, $zmin, $xmax, $ymax, $zmax);
					print OUT <<EOF;

(* Grey XY grid line, horizontal *)
newcurve linethickness .1 linetype solid marktype none gray 0.1
	pts $x1 $y1 $x2 $y2
EOF
				}

				# Since the grey grid was successful, let's try it on the
				# YZ plane.
				foreach my $yval (sort {$a <=> $b} keys %yvalues)
				{
					next if ($yval == $ymin);
					my ($x1, $y1)
						= &D3toD2($xmin, $yval, $zmin,
							$xmin, $ymin, $zmin, $xmax, $ymax, $zmax);
					my ($x2, $y2)
						= &D3toD2($xmin, $yval, $zmax,
							$xmin, $ymin, $zmin, $xmax, $ymax, $zmax);
					print OUT <<EOF;

(* Grey YZ grid line, vertical *)
newcurve linethickness .1 linetype solid marktype none gray 0.1
	pts $x1 $y1 $x2 $y2
EOF
				}
				foreach my $zval (sort {$a <=> $b} keys %zvalues)
				{
					next if ($zval == $zmin);
					my ($x1, $y1)
						= &D3toD2($xmin, $ymin, $zval,
							$xmin, $ymin, $zmin, $xmax, $ymax, $zmax);
					my ($x2, $y2)
						= &D3toD2($xmin, $ymax, $zval,
							$xmin, $ymin, $zmin, $xmax, $ymax, $zmax);
					print OUT <<EOF;

(* Grey YZ grid line, horizontal *)
newcurve linethickness .1 linetype solid marktype none gray 0.1
	pts $x1 $y1 $x2 $y2
EOF
				}

				# And finally, on the XZ plane, except that in this case,
				# we'll push it behind so it doesn't occlude real plots.
				foreach my $xval (sort {$a <=> $b} keys %xvalues)
				{
					next if ($xval == $xmin);
					my ($x1, $y1)
						= &D3toD2($xval, $ymax, $zmin,
							$xmin, $ymin, $zmin, $xmax, $ymax, $zmax);
					my ($x2, $y2)
						= &D3toD2($xval, $ymax, $zmax,
							$xmin, $ymin, $zmin, $xmax, $ymax, $zmax);
					print OUT <<EOF;

(* Grey XZ grid line, vertical *)
newcurve linethickness .1 linetype solid marktype none gray 0.1
	pts $x1 $y1 $x2 $y2
EOF
				}
				foreach my $zval (sort {$a <=> $b} keys %zvalues)
				{
					next if ($zval == $zmin);
					my ($x1, $y1)
						= &D3toD2($xmin, $ymax, $zval,
							$xmin, $ymin, $zmin, $xmax, $ymax, $zmax);
					my ($x2, $y2)
						= &D3toD2($xmax, $ymax, $zval,
							$xmin, $ymin, $zmin, $xmax, $ymax, $zmax);
					print OUT <<EOF;

(* Grey XZ grid line, horizontal *)
newcurve linethickness .1 linetype solid marktype none gray 0.1
	pts $x1 $y1 $x2 $y2
EOF
				}

				# Next, let's go through the points and print surfaces. The
				# surfaces will be printed as a grid, which means we have
				# to draw the "vertical" lines, i.e., lines connecting
				# points that have the same X value, and then we have to
				# draw the "horizontal" lines, i.e., lines connecting
				# points that have the same Y value. Also, for the
				# intermediate curves that form the surface, we shouldn't
				# do any labelling or else there'll be clutter.
				foreach my $plot (sort keys %Plots)
				{
					# Set the line, mark and colour for this curve set.
					my $linetype = $linetypes[$linetypeindex];
					$linetypeindex = ($linetypeindex + 1) % ($#linetypes + 1);
					my $marktype = $marktypes[$marktypeindex];
					$marktypeindex = ($marktypeindex + 1) % ($#marktypes + 1);
					my $colour = $colours[$colourindex];
					$colourindex = ($colourindex + 1) % ($#colours + 1);

					# Draw vertical lines. Also, prepare the data structure
					# for the horizontal lines.
					my %transpose = ();
					# We'll use the points below to generate the legend.
					my $labelx = undef; my $labely = undef;
					foreach my $xval (keys %{ $Plots{$plot} })
					{
						print OUT <<EOF;

(* Curve for plot $plot, vertical *)
newcurve linethickness $linethickness linetype $linetype
	marktype $marktype color $colour
	pts
EOF
						foreach my $yval (sort {$a <=> $b} keys %{ $Plots{$plot}{$xval} })
						{
							my $zval = $Plots{$plot}{$xval}{$yval};
							my ($xreal, $yreal)
								= &D3toD2($xval, $yval, $zval,
									$xmin, $ymin, $zmin, $xmax, $ymax, $zmax);
							print OUT "\t\t$xreal $yreal (* [$xval $yval $zval] *)\n";
							if (!defined($transpose{$yval}))
							{
								my %xhash = ();
								$transpose{$yval} = { %xhash };
							}
							$transpose{$yval}{$xval} = $zval;
						}
					}
					# Draw horizontal lines.
					foreach my $yval (keys %transpose)
					{
						print OUT <<EOF;

(* Curve for plot $plot, horizontal *)
newcurve linethickness $linethickness linetype $linetype
	marktype $marktype color $colour
	pts
EOF
						foreach my $xval (sort {$a <=> $b} keys %{ $transpose{$yval} })
						{
							my $zval = $transpose{$yval}{$xval};
							my ($xreal, $yreal)
								= &D3toD2($xval, $yval, $zval,
									$xmin, $ymin, $zmin, $xmax, $ymax, $zmax);
							print OUT "\t\t$xreal $yreal (* [$xval $yval $zval] *)\n";
							$labelx = $xreal;
							$labely = $yreal;
						}
					}
					# Here, we generate the label for this set of curves,
					# by drawing a dummy curve overlapping a genuine point.
					# It turns out that that point is the last point in the
					# set of curves simply because that's how we assigned
					# the labelx and labely points.
					print OUT <<EOF;

(* Legend for plot $plot *)
newcurve linethickness $linethickness linetype $linetype
	marktype $marktype color $colour label : $plot
	pts $labelx $labely
EOF
				}
			}

			close OUT;

			# Post-process if necessary to get the proper file.
			if ($realfile =~ /\.ps$/)
			{
				system("jgraph < $outfile > $realfile");
			}
			elsif ($realfile =~ /\.gif$/ or $realfile =~ /\.jpg$/)
			{
				# Make up the intermediate PS file name.
				my $tmpfile = $realfile;
				$tmpfile =~ s/\.gif$/.ps/;
				$tmpfile =~ s/\.jpg$/.ps/;
				system("jgraph < $outfile > $tmpfile");
				system("convert $tmpfile $realfile");
			}
		}
		else
		{
			die "Unsupported extension in output file $outfile.\n";
		}
	}
}