#!/usr/bin/perl -w

use strict;

unless ($ARGV[1]) {
  print "Usage: interp.pl exact_solution file1 [file2 [...]]\n";
  print "If exact==0, each series is just resampled to highest rate.\n";
  print "Otherwise, average relative error and deviation is output.\n";
  exit 1;
}

my $exact = $ARGV[0];

my $infinity = 1000000000000000000000000;

# data lines in each file should be in form:

# header stuff
# ...
# 1 distance [other fields...]
# n ...
# n ...
# (blank)

my (@n_series, @val_series);

# read in data series
for (my $argnum = 1; $argnum <= $#ARGV; ++$argnum) {
  my $ss = $argnum - 1;
  open IN, $ARGV[$argnum] or die "can't open ".$ARGV[$argnum].": $!\n";
  my $headerstuff = 1;
  my $i = 0;
  while (my $line = <IN>) {
	chomp($line);
	$line =~ s/^\s+//;
	$line =~ s/\s+$//;
	if (($headerstuff) and ($line =~ /^1/)) {$headerstuff = 0;}
	next if $headerstuff;
	if ($line =~ /^$/) {last;}
	my ($n, $val) = split(/\s+/, $line);
	$n = scalar($n);
	$n_series[$ss][$i] = $n;
	$val_series[$ss][$i] = $val;
	++$i;
  }
  close IN;
}

# diagnostics
#for (my $ss=0; $ss <= $#n_series; ++$ss) {
#  print "series $ss:\n";
#  for (my $i=0; $i < scalar @{$n_series[$ss]}; ++$i) {
#	print "$i\t" . $val_series[$ss][$i] . "\n";
#  }
#}

# interpolate and print
my (@i_series);
for (my $s=0; $s <= $#n_series; ++$s) {
  $i_series[$s] = 0;
}
my $all_done = 0;
data_loop: while (1) {
  # find next lowest n
  my $thisn = $infinity;
  for (my $s=0; $s <= $#n_series; ++$s) {
	if ($n_series[$s][$i_series[$s]] < $thisn) {
	  $thisn = $n_series[$s][$i_series[$s]];
	}
  }
  if ($thisn == $infinity) {exit;}
  # interpolate values at this n
  my (@interp);
  for (my $s=0; $s <= $#n_series; ++$s) {
	if ($n_series[$s][$i_series[$s]] == $thisn) {
	  # no need to interpolate
	  $interp[$s] = $val_series[$s][$i_series[$s]];
	} else {
	  # assume not at endpoint
	  my $lastn = $n_series[$s][$i_series[$s]-1];
	  my $nextn = $n_series[$s][$i_series[$s]];
	  my $lastval = $val_series[$s][$i_series[$s]-1];
	  my $nextval = $val_series[$s][$i_series[$s]];
	  my $grad = ($nextval-$lastval)/($nextn-$lastn);
	  $interp[$s] = $lastval + $grad * ($thisn-$lastn);
	}
  }
  if ($exact) {
	# calculate average relative error
	my $avg_val = 0;
	for (my $s=0; $s <= $#n_series; ++$s) {
	  $interp[$s] = 100 * ($interp[$s] - $exact) / $exact;
	  $avg_val += $interp[$s];
	}
	$avg_val /= $#n_series + 1;
	# calculate mean absolute deviation
	my $avg_dev = 0;
	for (my $s=0; $s <= $#n_series; ++$s) {
	  $avg_dev += abs($interp[$s] - $avg_val);
	}
	$avg_dev /= $#n_series + 1;
	# print this sample
	print "$thisn\t$avg_val\t$avg_dev\n";
  } else {
	# just print sample for each series
	print "$thisn";
	for (my $s=0; $s <= $#n_series; ++$s) {
	  print "\t".$interp[$s];
	}
	print "\n";
  }
  # increment indices of all series with this n
  for (my $s=0; $s <= $#n_series; ++$s) {
	if ($n_series[$s][$i_series[$s]] == $thisn) {
	  ++$i_series[$s];
	  if ($i_series[$s] == scalar(@{$n_series[$s]})) {
		# reached end of a data series
		last data_loop;
	  }
	}
  }
}