#!/usr/bin/perl -w
# renormalization group-based sampling error estimates:
# H. Flyvbjerg and H. G. Petersen, J. Chem. Phys. 91, 461 (1989)
# (C) Ilya Balabin, 2006

# average
sub average {
  my $s = 0.0;
  foreach $v (@_) {
    $s += $v;
  };
  return $s/($#_+1);
};

# standard deviation
sub stddev {
  my $ds = 0.0;
  my $s0 = average(@_);
  foreach $v (@_) {
    $ds += ($v-$s0)**2;
  };
  return sqrt($ds/($#_+1));
};

# sampling error estimates
sub serrs {
  my $se = stddev(@_) / sqrt($#_);
  my $dse = $se / sqrt(2 * $#_);
  return ($se, $dse);
};

# blocking transformation
sub block {
  my @vec = ();
  while ($#_>0) {
    my $v1 = shift;
    my $v2 = shift;
    push(@vec, 0.5*($v1+$v2));
  };
  return @vec;
};

# DEBUGGING: print vector elements in a single line
sub printvec {
  foreach $v (@_) {chomp $v; print "$v, ";};
  print "\n";
};

#######################################################
# usage syntax
$#ARGV >= 0 and ($v = shift) =~ m/^-h.*$/i ?
  die("Usage: renerr.pl file1.dat file2.dat ...\n") :
  unshift(@ARGV, $v);

# read input
@lines = <>;
die("Need at least 4 lines of single-column data\n")
  unless $#lines > 2;

# momenta
$x0 = average(@lines);
$dx0 = stddev(@lines);

# blocking transformations
$i = 0;
@serrs = ();
while ($#lines > 0) {
  push(@serrs, [ ($i, serrs(@lines)) ] );
  @lines = block(@lines);
  $i++;
};
$n = $i;

# print report
$x0 = sprintf("%10.5g", $x0);
$dx0 = sprintf("%10.5g", $dx0);
$line = "###############################################################\n";
print $line;
print "# Data average: $x0; Standard deviation: $dx0\n";
print $line;
print "# Blocking #\t\tsampling error\t\terror variation\n";
foreach $serr (@serrs) {
  foreach $v (@$serr) {
    printf("%-10.5g\t\t", $v);
  };
  print "\n";
};

# attempt to make an xmgrace plot
$xmgrace = `which xmgrace`;
chomp $xmgrace;
if (length($xmgrace)) {

  # scaling
  $ymin = $ymax = $serrs[0][1];
  foreach $serr (@serrs) {
    $ymin < $serr->[1] - $serr->[2] or $ymin = $serr->[1] - $serr->[2];
    $ymax > $serr->[1] + $serr->[2] or $ymax = $serr->[1] + $serr->[2];
  };

  # parameters
  $q1 = "'"; $q2 = '"';
  @pexecs = ();
  push(@pexecs, "${q1}SUBTITLE ${q2}\\z{1.1}Data average: ${x0}; Standard deviation: ${dx0}${q2}${q1}");
  push(@pexecs, "${q1}SUBTITLE FONT ${q2}Helvetica${q2}${q1}");
  push(@pexecs, "${q1}XAXIS LABEL ${q2}\\z{1.1}Blocking transformation #${q2}${q1}");
  push(@pexecs, "${q1}XAXIS LABEL FONT ${q2}Helvetica${q2}${q1}");
  push(@pexecs, "${q1}XAXIS TICKLABEL FONT ${q2}Helvetica${q2}${q1}");
  push(@pexecs, "${q1}YAXIS LABEL ${q2}\\z{1.1}Sampling error estimate${q2}${q1}");
  push(@pexecs, "${q1}YAXIS LABEL FONT ${q2}Helvetica${q2}${q1}");
  push(@pexecs, "${q1}YAXIS TICKLABEL FONT ${q2}Helvetica${q2}${q1}");
  $viewport = "-viewport 0.175 0.1 1.25 0.9";
  $world = "-world 0 $ymin $n $ymax";

  # command line arguments
  $xmstring = "| $xmgrace -pipe -noask -free $viewport $world -settype xydy -autoscale xy -settype xydy";
  foreach $p (@pexecs) {
    $xmstring = "$xmstring -pexec $p";
  };

  # pipe
  open FH, "$xmstring"
    or die("Can not pipe to ${xmgrace}: $! \n");

  # data
  foreach $serr (@serrs) {
    foreach $v (@$serr) {
      print FH "$v   ";
    };
    print FH "\n";
  };
  close FH or die("Can not close pipe to ${xmgrace}: $! \n");
} else {
  warn("WARNING: Xmgrace not found in the path, sampling errors not plotted.\n")
};

__END__