package Math;

#
# TODO
#
# 02  configurable analyze*() functions
#     - which confidence intervals and quantiles to calculate
#     - let analyze* call function to add confidence intervals and quantiles to stats
#     - separate call possible
# 05  write pod
# 06  remote deprecated functions
# 07  return gnuplot data instead of writing it to a file
#     (plot_histogram, plot_distribution)
# 08  write build_residual_from_values
# 09  write ensure_width()
# 10  use Carp
# 11  how else to call residual
#
#
# DONE
#
#   unknown
#     rename histo to histogram in functions and comments
#     revise visualize_residual_full() to use analyze_histo()
#     update @EXPORT with new names
#

### PERL INTERNAL

use strict;
use warnings;
use English;

### PUBLIC VARIABLES

our $VERSION = '1.0';

### INHERITANCE

require Exporter;
our @ISA = (
	'Exporter'
);

### EXPORTS

our @EXPORT = (
	'mergesort',
	'quicksort',
	
	'build_histogram',
	'expand_histogram',
	'normalize_histogram',
	'build_distribution',
	'build_distribution_smaller',
	'build_distribution_greater',
	'calculate_histogram_error',
	'calculate_distribution_error',
	'build_residual_from_histogram',

	'power',
	'round',

	'cut_above_quantile',
	'cut_below_quantile',

	'analyze',
	'analyze_sort',
	'analyze_histogram',
	'get_quantile',

	'analyze_basic',
	'analyze_quantile',
	'analyze_confidence',

	'precision',
	'display',
	'display_most',
	'plot_histogram_from_values',
	'plot_histogram',
	'plot_distribution_from_values',
	'plot_distribution',
	'plot_residual',
);

### PUBLIC METHODS

=pod

=head1 NAME

Math - statistical analysis

=head1 DESCRIPTION

this module provides functions for the statistical analysis of series of values.

=head1 PUBLIC METHODS

=over

=cut

# mergesort(array, start, end) {{{
=pod

=item mergesort(array, start, end)

sort I<array> inplace

=cut
sub mergesort {
	my $ref_data = shift;
	my $l = shift;
	my $r = shift;

	if ($r - $l > 0) {
		my $m = int(($l + $r) / 2);

		&mergesort($ref_data, $l, $m);
		&mergesort($ref_data, $m + 1, $r);

		my $i = $l;
		my $j = $m + 1;
		my $temp;
		while ($i <= $m and $j <= $r) {
			if ($ref_data->[$i] > $ref_data->[$j]) {
				$temp = $ref_data->[$j];
				for (my $k = $m; $k >= $i; --$k) {
					$ref_data->[$k + 1] = $ref_data->[$k];
				}
				++$m;
				$ref_data->[$i] = $temp;

				++$j;
			}
			++$i;
		}
	}
}
# }}}
# quicksort(array, start, end) {{{
=pod

=item quicksort(array, start, end)

sort I<array> inplace

=cut
sub quicksort {
	my $ref_data = shift;
	my $p = shift;
	my $r = shift;

	if ($p < $r) {
		my $temp;
		my $q;

		my $x = $ref_data->[$p];
		my $i = $p - 1;
		my $j = $r + 1;
		while (1) {
			do {
				--$j;
			} until ($ref_data->[$j] <= $x);
			do {
				++$i;
			} until ($ref_data->[$i] >= $x);
			if ($i < $j) {
				$temp = $ref_data->[$i];
				$ref_data->[$i] = $ref_data->[$j];
				$ref_data->[$j] = $temp;

			} else {
				$q = $j;
				last;
			}
		}

		&quicksort($ref_data, $p, $q);
		&quicksort($ref_data, $q + 1, $r);
	}
}
# }}}

# build_histogram(values, granularity, size): histogram {{{
=pod

=item build_histogram(values, granularity, size): histogram

build histogram from I<values> in I<size> buckets of I<granularity> width.

=cut
sub build_histogram {
	my $values = shift;
	my $gran = shift;
	my $size = shift;

	my $histo = [];
	for (my $i = 0; $i < $size; ++$i) {
		$histo->[$i] = 0;
	}

	foreach my $value (@$values) {
		my $i = int($value / $gran);
		if (not $i < $size) {
			$i = $size - 1;
		}
		$histo->[$i] += 1;
	}

	return $histo;
}
# }}}
# expand_histogram(histogram, granularity, offset, start, stop): values {{{
=pod

=item expand_histogram(histogram, granularity, offset, start, stop): values

generate values from I<histogram> assuming a bucket width of I<granularity> with initial bucket value I<offset> starting with bucket number I<start> up to bucket number I<stop>.

=cut
sub expand_histo {
	my $histo = shift;
	my $gran = shift;
	my $offset = shift;
	my $start = shift;
	my $stop = shift;

	if (not defined($offset)) {
		$offset = 0;
	}
	if (not defined($start)) {
		$start = 0;
	}
	if (not defined($stop)) {
		$stop = @$histo;
	}

	my $values = [];
	for (my $i = $start; $i < $stop; ++$i) {
		if ($histo->[$i] > 0) {
			for (my $j = 0; $j < $histo->[$i]; ++$j) {
				push @$values, $offset + ($i * $gran);
			}
		}
	}

	return $values;
}
# }}}
# normalize_histogram(histogram) {{{
=pod

=item normalize_histogram(histogram)

normalize I<histogram> with respect to the total number of values. this normalization is performed in place.

=cut
sub normalize_histogram {
	my $histo = shift;

	my $sum = 0;
	foreach my $bucket (@$histo) {
		$sum += $bucket;
	}

	for (my $i = 0; $i < @$histo; ++$i) {
		$histo->[$i] /= $sum;
	}
}
# }}}
# build_distribution(histogram): distribution {{{
=pod

=item build_distribution(histogram): distribution

alias for build_distribution_smaller

=cut
sub build_distribution {
	return &build_distribution_smaller(@_);
}
# }}}
# build_distribution_smaller(histogram): distribution {{{
=pod

=item build_distribution_smaller(histogram): distribution

build a distribution from a normalized I<histogram>.

=cut
sub build_distribution_smaller {
	my $histo = shift;

	&normalize_histogram($histo);

	my $distribution = [];
	my $sum = 0;
	foreach my $bucket (@$histo) {
		$sum += $bucket;
		push @$distribution, $sum;
	}

	return $distribution;
}
# }}}
# build_distribution_greater(histogram): distribution {{{
=pod

=item build_distribution_greater(histogram): distribution

build a distribution from a normalized I<histogram>.

=cut
sub build_distribution_greater {
	my $histo = shift;

	&normalize_histogram($histo);

	my $distribution = [];
	my $sum = 1;
	foreach my $bucket (@$histo) {
		$sum -= $bucket;
		push @$distribution, $sum;
	}

	return $distribution;
}
# }}}
# calculate_histogram_error(ref_histo, cur_histo): histogram {{{
=pod

=item calculate_histogram_error(ref_histo, cur_histo): histogram

calculate the per bucket error of histogram I<cur_histo> with respect to histogram I<ref_histo>.

=cut
sub calculate_histogram_error {
	my $ref_histo = shift;
	my $cur_histo = shift;

	my $buckets = @$ref_histo;
	my $diff_abs_histo = [];
	for (my $i = 0; $i < $buckets; ++$i) {
		$diff_abs_histo->[$i] = $cur_histo->[$i] - $ref_histo->[$i];
	}

	return $diff_abs_histo;
}
# }}}
# build_residual_from_histogram(histogram, granularity) {{{
=pod

=item residual_from_histogram(histogram, granularity)

TEXT

=cut
sub build_residual_from_histogram {
	my $histo = shift;
	my $gran = shift;

	my $residual = [];
	for (my $i = 0; $i < $#$histo; ++$i) {
		my $stats = &analyze_histogram($histo, $gran, $i);
		push @$residual, $stats;
	}

	return $residual;
}
# }}}

# power(value, exponent): value {{{
=pod

=item power(value, exponent): value 

calculate I<value> to the power of I<exponent>.

=cut
sub power {
	my $value = shift;
	my $exponent = shift;

	my $result = 1;
	for (my $i = 1; $i <= $exponent; ++$i) {
		$result *= $value;
	}

	return $result;
}
# }}}
# round(value, digits): value {{{
=pod

=item round(value, digits): value

mathematically round I<value> using I<digits> significant digits.

=cut
sub round {
	my $value = shift;
	my $digits = shift;

	if (not defined($digits)) {
		$digits = 0;
	}
	$value *= &power(10, $digits);

	my $rest = $value - int($value);
	if ($rest > 0.4) {
		$value = int($value) + 1;

	} else {
		$value = int($value);
	}

	return $value / &power(10, $digits);
}
# }}}

# cut_above_quantile(values, p): value {{{
=pod

=item cut_above_quantile(values, p): value

BLARG

=cut
sub cut_above_quantile {
	my $values = shift;
	my $p = shift;

	my $n = @$values;
	
	my $i = int($p * $n);
	while (@$values > $i) {
		pop @$values;
	}
}
# }}}
# cut_below_quantile(values, p): value {{{
=pod

=item cut_below_quantile(values, p): value

BLARG

=cut
sub cut_below_quantile {
	my $values = shift;
	my $p = shift;

	my $n = @$values;
	
	my $i = int($p * $n) + 1;
	while (@$values > $i) {
		shift @$values;
	}
}
# }}}

# analyze(values): stats{{{
=pod

=item analyze(values): stats

analyze I<values> returning the number, sum, minimum and maximum of values as well as the mean value, variance, and mean deviation. also the confidence interval for a 90%, 95% and 99% confidence and the 0.1, 0.25, 0.5, 0.75, 0.9 quantiles are added.

=cut
sub analyze {
	my $values = shift;

	my $stats = {
		'num' => 0,
		'sum' => 0,
	
		'min' => 0,
		'max' => 0,

		'avg' => 0,
		'var' => 0,
		'dev' => 0,

		'c90' => 0,
		'c95' => 0,
		'c99' => 0,

		'q10' => 0,
		'q25' => 0,
		'q50' => 0,
		'q75' => 0,
		'q90' => 0
	};

	if (not defined($values) or not @$values) {
		return $stats;
	}

	$stats->{'num'} = @$values;
	$stats->{'min'} = $values->[0];
	$stats->{'max'} = $values->[$#$values];

	for (my $i = 0; $i < @$values; ++$i) {
		my $value = $values->[$i];
		
		$stats->{'sum'} += $value;
	}
	$stats->{'avg'} = $stats->{'sum'} / $stats->{'num'};

	for (my $i = 0; $i < @$values; ++$i) {
		my $value = $values->[$i];

		$stats->{'var'} += &power($stats->{'avg'} - $value, 2);
	}
	$stats->{'var'} /= $stats->{'num'};
	$stats->{'dev'} = sqrt($stats->{'var'});

	# confidence interval
	#
	# K_s: Konfidenzinterval
	# x: Mittelwert
	# z: 
	#    90% = 1.645
	#    95% = 1.96
	#    99% = 2.58
	# s: Standardabweichung
	# n: Stichproben
	#
	# K_s = x +/- z s / sqrt(n)
	my $temp = $stats->{'dev'} / sqrt($stats->{'num'});
	$stats->{'c90'} = 1.645 * $temp;
	$stats->{'c95'} = 1.960 * $temp;
	$stats->{'c99'} = 2.580 * $temp;

	$stats->{'q10'} = &get_quantile($values, 0.10);
	$stats->{'q25'} = &get_quantile($values, 0.25);
	$stats->{'q50'} = &get_quantile($values, 0.50);
	$stats->{'q75'} = &get_quantile($values, 0.75);
	$stats->{'q90'} = &get_quantile($values, 0.90);

	return $stats;
}
# }}}
# analyze_sort(values): stats {{{
=pod

=item analyze_sort(values): stats

the list I<values> is sorted inplace, analyzed by I<analyze()> and the result returned.

=cut
sub analyze_sort {
	my $values = shift;

	if (defined($values)) {
		&quicksort($values, 0, $#$values);
	}
	return &analyze($values);
}
# }}}
# analyze_histogram(histogram, granularity, offset, start, stop): stats {{{
=pod

=item analyze_histogram(histogram, granularity, offset, start, stop): stats

this function analyzes a histogram with bucket width I<granularity> and initial bucket value I<offset> starting with bucket number I<start> up to bucket number I<stop>.

=cut
sub analyze_histogram {
	my $histo = shift;
	my $gran = shift;
	my $offset = shift;
	my $start = shift;
	my $stop = shift;

	if (not defined($offset)) {
		$offset = 0;
	}
	if (not defined($start)) {
		$start = 0;
	}
	if (not defined($stop)) {
		$stop = @$histo;
	}

	my $stats = {
		'num' => 0,
		'sum' => 0,
	
		'min' => 0,
		'max' => 0,

		'avg' => 0,
		'var' => 0,
		'dev' => 0,

		'c90' => 0,
		'c95' => 0,
		'c99' => 0,

		'q10' => 0,
		'q25' => 0,
		'q50' => 0,
		'q75' => 0,
		'q90' => 0
	};

	if (not defined($histo) or not @$histo) {
		return $stats;
	}

	$stats->{'min'} = $offset;
	$stats->{'max'} = $offset + @$histo * $gran;

	for (my $i = 0; $i < @$histo; ++$i) {
		my $value = $histo->[$i];
		
		$stats->{'num'} += $value;
		$stats->{'sum'} += $value * ($offset + $i * $gran);
	}
	$stats->{'avg'} = $stats->{'sum'} / $stats->{'num'};

	for (my $i = 0; $i < @$histo; ++$i) {
		my $value = $histo->[$i];

		$stats->{'var'} += $value * &power($stats->{'avg'} - ($offset + $i * $gran), 2);
	}
	$stats->{'var'} /= $stats->{'num'};
	$stats->{'dev'} = sqrt($stats->{'var'});

	# confidence interval
	#
	# K_s: Konfidenzinterval
	# x: Mittelwert
	# z: 
	#    90% = 1.645
	#    95% = 1.96
	#    99% = 2.58
	# s: Standardabweichung
	# n: Stichproben
	#
	# K_s = x +/- z s / sqrt(n)
	my $temp = $stats->{'dev'} / sqrt($stats->{'num'});
	$stats->{'c90'} = 1.645 * $temp;
	$stats->{'c95'} = 1.960 * $temp;
	$stats->{'c99'} = 2.580 * $temp;

	$stats->{'q10'} = &get_quantile_from_histogram($histo, $gran, $start, $stop, $stats->{'num'}, 0.10);
	$stats->{'q25'} = &get_quantile_from_histogram($histo, $gran, $start, $stop, $stats->{'num'}, 0.25);
	$stats->{'q50'} = &get_quantile_from_histogram($histo, $gran, $start, $stop, $stats->{'num'}, 0.50);
	$stats->{'q75'} = &get_quantile_from_histogram($histo, $gran, $start, $stop, $stats->{'num'}, 0.75);
	$stats->{'q90'} = &get_quantile_from_histogram($histo, $gran, $start, $stop, $stats->{'num'}, 0.90);

	return $stats;
}
# }}}
# get_quantile(values, p): p-quantile {{{
=pod

=item get_quantile(values, p): p-quantile

returns the I<p>-quantile from sorted I<values>

=cut
sub get_quantile {
	my $values = shift;
	my $p = shift;

	my $n = @$values;

	my $i = $p * $n;
	# the bucket numer $i is not integer: round up
	if ($i != int($i)) {
		return $values->[int($i)]; # rounding up is done implizitly because array indices run from 0 to n-1

	# the bucket number $i is an integer: mean value of bucket number $i and $i+1
	} else {
		return ($values->[$i - 1] + $values->[$i]) / 2; # remember: array indices run from 0 to n-1
	}
}
# }}}
# get_quantile_from_histogram(histogram, p): p-quantile {{{
=pod

=item get_quantile_from_histogram(histogram, p): p-quantile

returns the I<p>-quantile from a I<histogram>

=cut
sub get_quantile_from_histogram {
	my $histo = shift;
	my $gran = shift;
	my $start = shift;
	my $stop = shift;
	my $n = shift;
	my $p = shift;

	if ($p > 1) {
		die 'impossible quantile specified: p=' . $p . ' but must be 0<=p<=1' . "\n";
	}

	my $i = $p * $n;

	my $num = 0;
	my $bucket = $start - 1;
	while ($num < $i) {
		++$bucket;
		$num += $histo->[$bucket];
	}

	if ($i != int($i)) {
		return $bucket * $gran;

	} else {
		# last entry in bucket (interpolate with last bucket)
		if ($i == $num) {
			return (($bucket - 1) * $gran + $bucket * $gran) / 2;

		# first entry in bucket (interpolate with next bucket)
		} elsif ($i == ($num - $histo->[$bucket] + 1)) {
			return ($bucket * $gran + ($bucket + 1) * $gran) / 2;

		} else {
			return $bucket * $gran;
		}
	}
}
# }}}

# analyze_basic(values): stats {{{
=pod

=item analyze_basic(values): stats

TEXT

=cut
sub analyze_basic {
	my $values = shift;

	my $stats = {
		'num' => 0,
		'sum' => 0,
	
		'min' => 0,
		'max' => 0,

		'avg' => 0,
		'var' => 0,
		'dev' => 0,
	};

	if (not defined($values) or not @$values) {
		return $stats;
	}

	$stats->{'num'} = @$values;
	$stats->{'min'} = $values->[0];
	$stats->{'max'} = $values->[$#$values];

	for (my $i = 0; $i < @$values; ++$i) {
		my $value = $values->[$i];
		
		$stats->{'sum'} += $value;
	}
	$stats->{'avg'} = $stats->{'sum'} / $stats->{'num'};

	for (my $i = 0; $i < @$values; ++$i) {
		my $value = $values->[$i];

		$stats->{'var'} += &power($stats->{'avg'} - $value, 2);
	}
	$stats->{'var'} /= $stats->{'num'};
	$stats->{'dev'} = sqrt($stats->{'var'});

	return $stats;
}
# }}}
# analyze_confidence(values, confidence_definitions, stats): stats {{{
=pod

=item analyze_confidence(values, confidence_definitions, stats): stats

TEXT

=cut
sub analyze_confidence {
	my $values = shift;
	my $confidences = shift;
	my $stats = shift;

	if (not defined($stats)) {
		$stats = {};
	}

	if (not exists($stats->{'dev'}) or not exists($stats->{'num'})) {
		die 'you need to run analyze_basic() first' . "\n";
	}

	if (not defined($values) or not @$values) {
		return $stats;
	}

	my $zvalues = {
		'50' => 0.6745,
		'80' => 1.28,
		'90' => 1.645,
		'95' => 1.96,
		'96' => 2.05,
		'98' => 2.33,
		'99' => 2.58
	};

	my $temp = $stats->{'dev'} / sqrt($stats->{'num'});
	foreach my $confidence (@$confidences) {
		my ($tag) = ($confidence =~ m/^0\.(\d+)$/);
		$tag .= '0' x (2 - length($tag));
		if (not exists($zvalues->{$tag})) {
			die 'unknown confidence' . "\n";
		}
		$stats->{'c' . $tag} = $zvalues->{$tag} * $temp;
	}

	return $stats;
}
# }}}
# analyze_quantiles(values, quantile_definitions, stats): stats {{{
=pod

=item analyze_quantiles(values, quantile_definitions, stats): stats

TEXT

=cut
sub analyze_quantiles {
	my $values = shift;
	my $quantiles = shift;
	my $stats = shift;

	if (not defined($stats)) {
		$stats = {};
	}

	if (not defined($values) or not @$values) {
		return $stats;
	}

	foreach my $quantile (@$quantiles) {
		my ($tag) = ($quantile =~ m/^0\.(\d+)$/);
		$tag = '0' x (2 - length($tag));
		$stats->{'q' . $tag} = &get_quantile($values, $quantile);
	}

	return $stats;
}
# }}}

# precision(stats, digits): stats {{{
=pod

=item precision(stats, digits): stats

adjust precision of all values in I<stats> to have I<digits> significant digits

=cut
sub precision {
	my $stats = shift;
	my $digits = shift;

	foreach my $key (keys %$stats) {
		$stats->{$key} = &round($stats->{$key}, $digits);
	}

	return $stats;
}
# }}}
# display(stats): string {{{
=pod

=item display(stats): string

returns string containing all common stats

=cut
sub display {
	my $stats = shift;

	if (keys %$stats == 0) {
		return 'EMPTY';
	}

	my $has_stats = 0;
	foreach my $key (keys %$stats) {
		if (defined($stats->{$key}) and $stats->{$key} != 0) {
			$has_stats = 1;
		}
	}
	if (not $has_stats) {
		return 'EMPTY';
	}

	#return 'num=' . $stats->{'num'} . ' ' . &display_most($stats);
	return $stats->{'num'} . ' ' . &display_most($stats);
}
# }}}
# display_most(stats): string {{{
=pod

=item display_most(stats): string

returns string containing most common stats

=cut
sub display_most {
	my $stats = shift;

	if (keys %$stats == 0) {
		return 'EMPTY';
	}

	my $has_stats = 0;
	foreach my $key (keys %$stats) {
		if (defined($stats->{$key}) and $stats->{$key} != 0) {
			$has_stats = 1;
		}
	}
	if (not $has_stats) {
		return 'EMPTY';
	}

	#return 'min=' . $stats->{'min'} . '/avg=' . $stats->{'avg'} . '/max=' . $stats->{'max'} . ' (dev=' . $stats->{'dev'} . '/c95=' . $stats->{'c95'} . ') [q10=' . $stats->{'q10'} . '/q25=' . $stats->{'q25'} . '/median=' . $stats->{'q50'} . '/q75=' . $stats->{'q75'} . '/q90=' . $stats->{'q90'} . ']';
	return $stats->{'min'} . '/' . $stats->{'avg'} . '/' . $stats->{'max'} . ' (' . $stats->{'dev'} . '/' . $stats->{'c95'} . ') [' . $stats->{'q10'} . '/' . $stats->{'q25'} . '/' . $stats->{'q50'} . '/' . $stats->{'q75'} . '/' . $stats->{'q90'} . ']';
}
# }}}
# plot_histogram_from_values(values, granularity, size, name) {{{
=pod

=item plot_histogram_from_values(values, granularity, size, name)

TEXT

=cut
sub plot_histogram_from_values {
	my $values = shift;
	my $gran = shift;
	my $size = shift;
	my $name = shift;

	if (not defined($gran) or not defined($size)) {
		die 'unable to plot histo' . "\n";
	}
	
	my $histo = &build_histogram($values, $gran, $size);
	&plot_histogram($histo, $gran, $name);
}
# }}}
# plot_histogram(histogram, granularity, name) {{{
=pod

=item plot_histogram(histogram, granularity, name)

TEXT

=cut
sub plot_histogram {
	my $histo = shift;
	my $gran = shift;
	my $name = shift;

	if (not defined($gran)) {
		die 'unable to plot histo' . "\n";
	}

	my $size = @$histo;
	
	my $handle;
	if (not open($handle, '>' . $name . '.gnuplot')) {
		die 'unable to open file <' . $name . '.gnuplot>' . "\n";
	}

	print $handle 'set terminal png small' . "\n";
	print $handle 'set output "' . $name . '.png"' . "\n";
	print $handle 'set title "histo (gran=' . $gran . ' / size=' . $size . ')"' . "\n";
	print $handle 'set autoscale' . "\n";
	print $handle 'set xlabel "bucket [s]"' . "\n";
	print $handle 'set xrange [0:]' . "\n";
	print $handle 'set ylabel "# of values"' . "\n";
	print $handle 'set yrange [0:]' . "\n";
	print $handle 'set ytics nomirror' . "\n";
	print $handle "\n";
	print $handle 'plot \\' . "\n";
	print $handle '\'-\' using 1:2 \\' . "\n";
	print $handle '     axes x1y1 \\' . "\n";
	print $handle '     notitle \\' . "\n";
	print $handle '     with boxes \\' . "\n";
	print $handle '          linetype 1 \\' . "\n";
	print $handle "\n";

	for (my $i = 0; $i < $size; ++$i) {
		print $handle '' . ($i * $gran + $gran / 2) . ' ' . $histo->[$i] . "\n";
	}
	print $handle 'e' . "\n";

	close($handle);
}
# }}}
# plot_distribution(distribution, granularity, name) {{{
=pod

=item plot_distribution(distribution, granularity, name)

TEXT

=cut
sub plot_distribution {
	my $distribution = shift;
	my $gran = shift;
	my $name = shift;

	if (not defined($gran)) {
		die 'unable to plot histo' . "\n";
	}
	$name .= '.gnuplot';

	my $size = @$distribution;

	my $handle;
	if (not open($handle, '>' . $name)) {
		die 'unable to open file <' . $name . '>' . "\n";
	}

	print $handle 'set terminal png small' . "\n";
	#print $handle 'set output "' . $basename . '.datastats.png"' . "\n";
	print $handle 'set title "histo (gran=' . $gran . ' / size=' . $size . ')' . "\n";
	print $handle 'set autoscale' . "\n";
	print $handle 'set xlabel "bucket [s]"' . "\n";
	print $handle 'set xrange [0:]' . "\n";
	print $handle 'set ylabel "probability"' . "\n";
	print $handle 'set yrange [0:]' . "\n";
	print $handle 'set ytics nomirror' . "\n";
	print $handle "\n";
	print $handle 'plot \\' . "\n";
	print $handle '\'-\' using 1:2 \\' . "\n";
	print $handle '     axes x1y1 \\' . "\n";
	print $handle '     notitle \\' . "\n";
	print $handle '     with linespoints \\' . "\n";
	print $handle '          linetype 1 \\' . "\n";
	print $handle "\n";

	for (my $i = 0; $i < $size; ++$i) {
		print $handle '' . $i . ' ' . $distribution->[$i] . "\n";
	}
	print $handle 'e' . "\n";

	close($handle);
}
# }}}
# plot_residual(residual, granularity, name) {{{
=pod

=item plot_residual(resitual, granularity, name)

TEXT

=cut
sub plot_residual {
	my $residual = shift;
	my $gran = shift;
	my $name = shift;
	
	my $handle;
	if (not open($handle, '>' . $name . '.gnuplot')) {
		die 'unable to open file <' . $name . '.gnuplot>' . "\n";
	}

	print $handle 'set terminal png small' . "\n";
	print $handle 'set output "' . $name . '.png"' . "\n";
	print $handle 'set title "residual of histo (gran=' . $gran . ')"' . "\n";
	print $handle 'set autoscale' . "\n";
	print $handle 'set xlabel "bucket [s]"' . "\n";
	print $handle 'set xrange [0:]' . "\n";
	print $handle 'set ylabel "residual"' . "\n";
	print $handle 'set yrange [0:]' . "\n";
	print $handle 'set ytics nomirror' . "\n";
	print $handle "\n";
	print $handle 'plot \\' . "\n";
	print $handle '\'-\' using 1:2 \\' . "\n";
	print $handle '     axes x1y1 \\' . "\n";
	print $handle '     notitle \\' . "\n";
	print $handle '     with lines \\' . "\n";
	print $handle '          linetype 3 \\' . "\n";
	print $handle '     , \\' . "\n";
	print $handle '\'-\' using 1:2:3 \\' . "\n";
	print $handle '     axes x1y1 \\' . "\n";
	print $handle '     notitle \\' . "\n";
	print $handle '     with yerrorbars \\' . "\n";
	print $handle '          linetype 3 \\' . "\n";
	print $handle '          pointtype 0 \\' . "\n";
	print $handle "\n";
	
	my @data = ();
	for (my $i = 0; $i < @$residual; ++$i) {
		my $stats = $residual->[$i];
		push @data, $i . ' ' . ($stats->{'avg'} - $i * $gran) . ' ' . $stats->{'c95'};
	}

	foreach my $entry (@data) {
		print $handle $entry . "\n";
	}
	print $handle 'e' . "\n";
	foreach my $entry (@data) {
		print $handle $entry . "\n";
	}
	print $handle 'e' . "\n";

	close($handle);
}
# }}}

=pod

=back

=head1 AUTHOR

Nicholas Dille E<lt>webmaster@rakshas.deE<gt>

=cut

1;