
$filename="";

do {
	show_summary();
	
	if ($filename eq "") {
		$command="load";
	}

	print "\nCommands {load, q, fltr, hist}\n";
	$command=read_input("What next?","$command");
	$command =~ tr/[A-Z]/[a-z]/;
	if ($command eq "load") {
		load();
	} elsif ($command eq "fltr") {
		set_filter();
	} elsif ($command eq "hist") {
		histogram();
	} 
} until ($command eq "q");

print "\nBye!\n\n\n";
die;

sub set_filter() {
	if ($isize <=0) {
		print "don't know how many columns there are\n";
		return;
	}
	
	$col_filter=read_input("Which column?",0);
	if (($col_filter >= $isize) || ($col_filter < 0)){
		print "$col_filter doesn't sound right\n";
		return;
	}
	
	print "Column $col_filter, min $min[$col_filter], max $max[$col_filter]\n";
	
	$filter_start=read_input("Enter start of range (inclusive of point)",$min[$col_filter]);
	$filter_end=read_input("Enter end of range (exclusive of point)",$max[$col_filter]);
	if ($filter_end < $filter_start) {
		print "the end shouldn't be greater than the start\n";
		return;
	}
	$filter_flag="on";
	
	return;
}

sub load() {
	my ($a, @b, $i, $old_filename, $start);
	my ($header);  
	local (*fp1);
	
	$start=1;
	
	$old_filename=$filename;
	($filename eq "") && ($filename = "../jul06/day11/run1/posn_rel.1");
	$filename=read_input("what file?",$filename);
	
	($filename ne $old_filename) && ($filter_flag = "off");
	if ($filename =~/ice_depth/) {
		$ignore=read_input("Doing ice depth ,ignore depth '99' or '9999'","yes");
	}

	open(fp1,$filename) || die;
	$a=<fp1>;
	print "first line\n";
	print $a;
	chomp($a);

	$mode="tab";
	@b=split/	/,$a;
	$isize=$#b+1;

	if ($isize<=1) {
		$mode="spc";
		@b=split/\s+/,$a;
		$isize=$#b+1;
	}
	if ($isize<1) {
		print "number of columns < 1, surely that isn't right?!\n";
		die;
	}

# test if value is text by testing if ==0 while not containing '0'
	$i=-1;
	do {
		++$i;
	} while (($b[$i] == 0) && !($b[$i] =~/0/) && ($i<$isize));
	
	if ($i==$isize) {
		$header=1;
		for ($i=0;$i<$isize;++$i) {
			$label[$i]=$b[$i];
		}
		$a=<fp1>;
		print "second line\n";
		print $a;
		chomp($a);
	} else {
		$header=0;
		for ($i=0;$i<$isize;++$i) {
			$label[$i]=$i;
		}
	}

	close(fp1);
	print "$isize columns, mode $mode\n";

	for ($i=0;$i<$isize;$i++) {
		$sum[$i]=0;
		$sum2[$i]=0;
	}

	$count=0;
	$line_count=0;

	open(fp1,$filename);
	for ($i=0;$i<$header;++$i) {
		$a=<fp1>;
	}
	while ($a=<fp1>) {
		(($line_count%10000)==0) && print ".";
		if ($mode eq "tab") {
			@b=split/	/,$a;
		} elsif ($mode eq "spc") {
			@b=split/\s+/,$a;
		}
		
		if (($filter_flag ne "on")  || 
		  (($filter_start <= $b[$col_filter]) && ($b[$col_filter] < $filter_end))) {
			if (($ignore ne "yes") || (($b[2]!=99) && ($b[2]!=9999))) {
				if ($start==1) {
					for ($i=0;$i<$isize;$i++) {
						$min[$i]=1.0*$b[$i];
						$max[$i]=1.0*$b[$i];
					}
					$start=0;
				}
				for ($i=0;$i<$isize;$i++) {
					$sum[$i]+=$b[$i];
					$sum2[$i]+=$b[$i]*$b[$i];
					if ($b[$i]<$min[$i]) {
						$min[$i]=1.0*$b[$i];
					}
					if ($b[$i]>$max[$i]) {
						$max[$i]=1.0*$b[$i];
					}
				}
				++$count;
			}
		}
		++$line_count;
	}
	print "\n";
	close(fp1);
}

sub histogram() {
	my($i, $ilat, $count, $a, @b);
	local (*fp1, *fp);
	my($occurences, $zsuggest, $zstd, $zct, $weight, $label);
	my($hist_min, $hist_step, $hist_size, $hist_max, $hist_weight, @hist_bin);
	
	if ($isize <=0) {
		print "Is a file loaded yet?\n";
		return;
	}
	
	print "Doing Histograms\n";
	print "Filter $filter_flag\n";
	print "File $filename ($mode)\n";
	$ihist=read_input("Which columns you would like to histogram? (-1 to exit)",-1);
	($ihist < 0) && return;
	($ihist >= $isize) && return;

	$hist_min=read_input("What is the minimum of range?",$min[$ihist]);

	$z100=($max[$ihist]-$min[$ihist])/100.0;
	$zsuggest=$z100;
	$zstd=$std_dev[$ihist]/2.5;
	$zstd<$zsuggest && ($zsuggest=$zstd);
	$zct=100.0/0.3/$line_count*$std_dev[$ihist];
	$zct<$zsuggest && ($zsuggest=$zct);

	$label=sprintf "What is the step size? (100 %g, s.d. %g, ct %g)",$z100, $zstd, $zct;
	$hist_step=read_input("$label",$zsuggest);
	$hist_step==0 && return;
	$zsuggest=($max[$ihist]-$min[$ihist])/$hist_step;

	$hist_size=read_input("What is the number of steps?",$zsuggest);
	$hist_max=$hist_min+$hist_step*$hist_size;

	$hist_weight=read_input("Weighting? (area/uniform)","uniform");
	if ($hist_weight eq "area") {
		print "Will work in relative area on surface\n";
		$ilat=read_input("Which column has the latitude (in degrees)?",-1);
		$ilat < 0 && return;
		$ilat >=$isize && return;
		$min[$ilat] < -90 && return;
		$max[$ilat] > 90 && return;
	}

	for ($i=0;$i<$hist_size;++$i) {
		$hist_bin[$i]=0;
	}

	open(fp1,$filename);
	$count=0;
	$line_count=0;
	while ($a=<fp1>) {
		(($line_count%10000)==0) && print ".";
	#for ($j=0;$j<10;$j++) {
		if ($mode eq "tab") {
			@b=split/	/,$a;
		} elsif ($mode eq "spc") {
			@b=split/\s+/,$a;
		}
		if ($hist_weight eq "uniform") {
			$weight=1;
		} elsif ($hist_weight eq "area") { 
			$weight=1.0*cos($b[$ilat]*3.14159265/180.0);
		} else {
			die;
		}

		if (($filter_flag ne "on")  || 
		  (($filter_start <= $b[$col_filter]) && ($b[$col_filter] < $filter_end))) {
			if (($ignore ne "yes") || (($b[2]!=99) && ($b[2]!=9999))) {
				if (($b[$ihist] < $hist_max) && ($b[$ihist] > $hist_min)) {
					$hist_bin[($b[$ihist]-$hist_min)/$hist_step]+=$weight;
					$count+=$weight;
				}
			}
		}
		++$line_count;
	}
	close(fp1);
	print "\n";

	if ($hist_weight eq "uniform") {
		$label="Occurences";
	} elsif ($hist_weight eq "area") { 
		$label="Rel.Area";
	} else {
		die;
	}

	open(fp,">stat_hist.txt") || die;
	print "Histogram of file $filename, column $ihist\n";
	print "Min $min[$ihist],	Max $max[$ihist],	Av $av[$ihist],	Std Dev $std_dev[$ihist]\n";
	print "Number of values in histogram $count\n";
	print "Bin	Mid_bin value	$label\n";
	print fp "Histogram of file $filename, column $ihist\n";
	print fp "Min $min[$ihist],	Max $max[$ihist], Av $av[$ihist],	Std Dev $std_dev[$ihist]\n";
	print fp "Number of values in histogram $count\n";
	print fp "Bin	Mid_bin value	$label\n";

	for ($i=0;$i<$hist_size;++$i) {
		$value=$hist_min+(0.5+$i)*$hist_step;
		print "$i	$value	$hist_bin[$i]\n";
		print fp "$i	$value	$hist_bin[$i]\n";
	}
	
	print "Total number of lines $line_count\n";
	print fp "Total number of lines $line_count\n";
	if ($hist_weight eq "area") {
		print "Area counted	$count\n";
		print fp "Area counted	$count\n";
	} else {
		print "Lines counted $count\n";
		print fp "Lines counted $count\n";
	}
	
	close(fp);

	return;
}

sub floor {
    my($number) = shift;
    if ($number == int($number)) {
	    return $number;
    } else {
    	return int($number - 1*($number < 0));
    }
}

sub read_input($$) {
	my ($a, $b, $c);
	
	$a=$_[0];
	$b=$_[1];
	
	print "$a ($b) ";
	$c=<>;
	
	chop($c);
	if ($c eq "") {
		$c=$b;
	}
	
	return $c;
}

sub show_summary() {
	my ($name, $iname);
	print "\nSummary $filename ($mode)\n";
	print "Total lines $line_count, lines included $count\n";
	print "Filter	$filter_flag\n";
	if ($filter_flag eq "on") {
		print "Filter:	column $col_filter,	start $filter_start, end $filter_end\n";
	}

	print "\n";
	if ($isize < 5) {
	#print results along row
		print "min\n";
		for ($i=0;$i<$isize;$i++) {
			print "$min[$i]	";
		}
		print "\n";
		print "max\n";
		for ($i=0;$i<$isize;$i++) {
			print "$max[$i]	";
		}
		print "\n";
		print "average\n";
		for ($i=0;$i<$isize;$i++) {
			$av[$i]=$sum[$i]/$count;
			print sprintf "%g	",$av[$i];
		}
		print "\n";
		print "standard deviation\n";
		($count==1) && die;
		for ($i=0;$i<$isize;$i++) {
			$std_dev[$i]=sqrt(($sum2[$i]-$sum[$i]*$sum[$i]/$count)/($count-1));
			print sprintf "%g	",$std_dev[$i];
		}
		print "\n";
	} else {
	#too many columns for single row, instead list down
		print "Col#	min	max	average	std_dev\n";
		for ($i=0;$i< $isize;++$i) {
			$av[$i]=$sum[$i]/$count;
			$std_dev[$i]=sqrt(($sum2[$i]-$sum[$i]*$sum[$i]/$count)/($count-1));
			$_=$label[$i];
			$iname=tr/a-z//;
			if ($iname > 5) {
				$name=substr($_,0,6);
			} else {
				$name=$_;
			}
			print sprintf "$name	$min[$i]	$max[$i]	%g	%g\n", $av[$i],$std_dev[$i];
		}
	}

	return;
}