#!/usr/bin/perl

$| = 1;	# autoflush

# %lines is asociative array (hash)
my %lines;
my $i, $j;
my $value = "";
my $output = "";

# add zeros to matrix
for ($i = 0; $i < 256; $i++) {
	for ($j = 0; $j < 256; $j++) {
		$value = "$i $j";
		if (!exists $lines{"$value"}) {
			$lines{"$value"} = 0;
		}
	}
}

while (<STDIN>) {

	chop;
	if (exists $lines{"$_"}) {
		$lines{"$_"} = $lines{"$_"} + 1;
	}
	else {
		$lines{"$_"} = 1;
	}
}

# print header of the input file for "R" language interpreter
print "#!/usr/bin/R BATCH\n";
print "\n";
print "# Use R (\"GNU S\") interpreter to produce graphs\n";
print "# http://www.r-project.org/\n";
print "\n";
print "# inicialization: create vector 'z' from experimental data\n";
print "z <- c(";

# print count for each entry in hash
#foreach $key (sort keys %lines) {
#	print $key, " ", $lines{"$key"}, "\n";
#}

## print output
for ($i = 4; $i < 156; $i++) {
	for ($j = 4; $j < 209; $j++) {
		$value = "$i $j";
		print $lines{"$value"}, ", "; 
	}
#	print "\n";
}

print ")\n";

print "cat(\"Number of points:\", length(z), \"\\n\")\n";
print "\n";
print "mat <- matrix(z, nrow=152,ncol=205, byrow=\"T\")\n";
print "\n";
print "x <- rep(4:155, length.out=dim(mat)[1])\n";
print "y <- rep(4:208, length.out=dim(mat)[2])\n";
print "\n";
print "cat(\"Plotting histogram3D-notfitted.eps\\n\")\n";
print "postscript(file = \"histogram3D-notfitted.eps\", horizontal = FALSE, onefile = FALSE)\n";
print "persp(mat, xlab=\"napatie U[V]\", ylab=\"prud I[A]\", zlab=\"pocetnost\", theta=-55, phi=20, d=0.8)\n";
print "dev.off()\n";
print "\n";
print "\n";
print "# Fitting & smoothing\n";
print "for (k in 1:5) {\n";
print "	cat(\"Processing\", k, \"\\n\")\n";
print "	# fitting over y-z plane\n";
print "	for (i in 1:length(x)) {\n";
print "		for (j in 1:(length(y) - 2)) {\n";
print "			if (mat[i,j+1] <= mat[i,j]) {\n";
print "				if (mat[i,j+1] <= mat[i,j+2]) {\n";
print "					mat[i,j+1] <- (mat[i,j] + mat[i,j+2] )/2\n";
print "				}\n";
print "			}\n";
print "		}\n";
print "	}\n";
print "\n";
print "	# fitting over x-z plane\n";
print "	for (j in 1:length(y)) {\n";
print "		for (i in 1:(length(x) - 2)) {\n";
print "			if (mat[i+1,j] <= mat[i,j]) {\n";
print "				if (mat[i+1,j] <= mat[i+2,j]) {\n";
print "					mat[i+1,j] <- (mat[i,j] + mat[i+2,j] )/2\n";
print "				}\n";
print "			}\n";
print "		}\n";
print "	}\n";
print "}\n";
print "for (i in 1:length(x)) {\n";
print "	for (j in 1:(length(y) - 2)) {\n";
print "				mat[i,j+1] <- (mat[i,j] + mat[i,j+2] )/2\n";
print "	}\n";
print "}\n";
print "for (j in 1:length(y)) {\n";
print "	for (i in 1:(length(x) - 2)) {\n";
print "				mat[i+1,j] <- (mat[i,j] + mat[i+2,j] )/2\n";
print "	}\n";
print "}\n";
print "\n";
print "# matrix is fitted, create graph...\n";
print "cat(\"Plotting histogram3D.eps\\n\")\n";
print "postscript(file = \"histogram3D.eps\", horizontal = FALSE, onefile = FALSE)\n";
print "persp(mat, xlab=\"napatie U[V]\", ylab=\"prud I[A]\", zlab=\"pocetnost\", theta=-55, phi=20, d=0.8)\n";
print "dev.off()\n";
print "\n";
print "#X11()\n";
print "#persp(mat, xlab=\"napatie U[V]\", ylab=\"prud I[A]\", zlab=\"pocetnost\", theta=-55, phi=20, d=0.8)\n";
print "\n";
print "\n";
print "# create new matrix with smallest range\n";
print "new <- matrix(NA, nrow=150, ncol=200)\n";
print "for (i in seq(1, 150, by=1)) {  for (j in seq(1, 200, by=1 )) { new[i,j] <-  mat[i,j] }}\n";
print "\n";
print "# 1/2 resolution\n";
print "cat(\"Plotting histogram3D-2.eps\\n\")\n";
print "small <- matrix(NA, nrow=75, ncol=100)\n";
print "for (i in seq(1, 75, by=1)) {  for (j in seq(1, 100, by=1 )) { small[i,j] <-  new[2*i,2*j] }}\n";
print "postscript(file = \"histogram3D-2.eps\", horizontal = FALSE, onefile = FALSE)\n";
print "persp(small, xlab=\"napatie U[V]\", ylab=\"prud I[A]\", zlab=\"pocetnost\", theta=-55, phi=20, d=0.8)\n";
print "dev.off()\n";
print "rm(small)\n";
print "\n";
print "# 1/3 resolution\n";
print "cat(\"Plotting histogram3D-3.eps\\n\")\n";
print "small <- matrix(NA, nrow=50, ncol=66)\n";
print "for (i in seq(1, 50, by=1)) {  for (j in seq(1, 66, by=1 )) { small[i,j] <-  new[3*i,3*j] }}\n";
print "postscript(file = \"histogram3D-3.eps\", horizontal = FALSE, onefile = FALSE)\n";
print "persp(small, xlab=\"napatie U[V]\", ylab=\"prud I[A]\", zlab=\"pocetnost\", theta=-55, phi=20, d=0.8)\n";
print "dev.off()\n";
print "rm(small)\n";
print "\n";
print "# 1/4 resolution\n";
print "cat(\"Plotting histogram3D-4.eps\\n\")\n";
print "small <- matrix(NA, nrow=37, ncol=50)\n";
print "for (i in seq(1, 37, by=1)) {  for (j in seq(1, 50, by=1 )) { small[i,j] <-  new[4*i,4*j] }}\n";
print "postscript(file = \"histogram3D-4.eps\", horizontal = FALSE, onefile = FALSE)\n";
print "persp(small, xlab=\"napatie U[V]\", ylab=\"prud I[A]\", zlab=\"pocetnost\", theta=-55, phi=20, d=0.8)\n";
print "dev.off()\n";
print "rm(small)\n";
print "\n";
print "# 1/5 resolution\n";
print "cat(\"Plotting histogram3D-5.eps\\n\")\n";
print "small <- matrix(NA, nrow=30, ncol=40)\n";
print "for (i in seq(1, 30, by=1)) {  for (j in seq(1, 40, by=1 )) { small[i,j] <-  new[5*i,5*j] }}\n";
print "postscript(file = \"histogram3D-5.eps\", horizontal = FALSE, onefile = FALSE)\n";
print "persp(small, xlab=\"napatie U[V]\", ylab=\"prud I[A]\", zlab=\"pocetnost\", theta=-55, phi=20, d=0.8)\n";
print "dev.off()\n";
print "rm(small)\n";
print "\n";
print "# contoured picture at display\n";
print "xtmp <- rep(4:155, length.out=152)\n";
print "ytmp <- rep(4:155, length.out=152)\n";
print "for (i in 1:152) { xtmp[i] <- 67.1 * i /256}\n";
print "for (i in 1:152) { ytmp[i] <- 529 * i /256}\n";
print "\n";
print "haluz <- matrix(NA, nrow=152, ncol=152)\n";
print "for (i in 1:152) {\n";
print "	for (j in 1:152) {\n";
print "		haluz[i,j] <- mat[i,j]\n";
print "	}\n";
print "}\n";
print "\n";
print "# contoured picture in EPS grafic\n";
print "cat(\"Plotting histogram-contour.eps\\n\")\n";
print "postscript(file = \"histogram-contour.eps\", horizontal = FALSE, onefile = FALSE)\n";
print "contour(xtmp,ytmp,haluz)\n";
print "title(xlab=\"napatie U[V]\", ylab=\"prud I[A]\", font=3)\n";
print "dev.off()\n";
print "\n";