1#!/usr/bin/env perl 2use POSIX qw/floor/; 3use strict; 4use FileHandle; 5use Getopt::Long; 6 7my %config = do "/gprojects/qmcpack/qmcpack/utils/setup-qmc-conf.pl"; 8my $gnuplot = $config{gnuplot}; 9 10my $epsfile; 11GetOptions('eps=s' => \$epsfile); 12 13my $file = shift; 14my $start = shift; 15(defined $start) || ($start = 0); 16my $end = shift @ARGV; 17if (!(defined $end)) { 18 $end = -1; 19} 20 21open(FILE, "$file"); 22my @filedata; 23while (defined (my $line = <FILE>)) { 24 my @tmp = split('\s+', $line); 25 push(@filedata, [@tmp]); 26} 27close(FILE); 28 29 30my $min; 31my $max; 32($min, $max) = getMinMax(\@filedata, 2, $start, $end); 33my $enTicString = getTicString($min,$max); 34 35my $popmin; 36my $popmax; 37($popmin, $popmax) = getMinMax(\@filedata, 5, $start, $end); 38my $popTicString = getTicString($popmin,$popmax); 39 40my $plotstring = "set title \"$file\"\n"; 41if ($epsfile) { 42 $plotstring .= "set term post color enhanced 20\n set output \"$epsfile\"\n"; 43} 44$plotstring .= "set multiplot\n set ylabel \"Energy (Hartrees)\"\n "; 45$plotstring .= "set origin 0,0.5\n set size 0.82,0.5\n"; 46$plotstring .= "set format x \"\"\n set ytics $enTicString \n unset y2tics\n set lmargin 12\n"; 47if ($end < 0) { 48 $plotstring .= "plot [$start:][$min:$max] \"$file\" u 1:2 w l notitle, \"\" u 1:5 lt 3 notitle w l \n"; 49} else { 50 $plotstring .= "plot [$start:$end][$min:$max] \"$file\" u 1:2 w l notitle, \"\" u 1:5 lt 3 notitle w l \n"; 51} 52$plotstring .= "set format x \n unset title\n set xlabel \"DMC step\"\n set ylabel \"Population (Walkers)\"\n"; 53$plotstring .= "set xtics\n set ytics $popTicString\n set bmargin 4\n"; 54if ($end < 0) { 55 $plotstring .= "set origin 0,0\n plot [$start:][$popmin:$popmax] \"\" u 1:5 w l notitle\n"; 56} else { 57 $plotstring .= "set origin 0,0\n plot [$start:$end][$popmin:$popmax] \"\" u 1:5 w l notitle\n"; 58} 59$plotstring .= "set format x \"\"\n"; 60$plotstring .= "set origin 0.8,0.5\n set size 0.2, 0.5\n set title \"histogram\"\n"; 61$plotstring .= "set lmargin 0\n unset xlabel\n unset ylabel \n set bmargin 1\n"; 62$plotstring .= "unset ytics \n set y2tics $enTicString \n set format y2 \"\"\n "; 63$plotstring .= getGPLHistString(\@filedata, 2, $start, $end); 64$plotstring .= "set origin 0.8,0.0\n set size 0.2, 0.5\n unset title \n set y2tics $popTicString \n"; 65$plotstring .= "set format y2 \"\"\n set xlabel \"counts\"\n set bmargin 4\n"; 66$plotstring .= getGPLHistString(\@filedata, 5, $start, $end); 67$plotstring .= "unset multiplot\n"; 68unless($epsfile) { 69 $plotstring .= "pause -1\n"; 70} 71 72open(GPL, "|$gnuplot"); 73GPL->autoflush(1); 74print GPL $plotstring; 75unless($epsfile) { 76 my $blah = <>; # Hack to leave graph up until user hits a key 77} 78close(GPL); 79 80 81############################################################################# 82 83sub getTicString { 84 use POSIX qw/floor ceil/; 85 my $low = shift; 86 my $high = shift; 87# print " $low $high\n"; 88 my $range = $high - $low; 89 my $scale = 10 ** floor(log($range)/log(10)); 90 my $ntics = int($range/$scale); 91 if ($ntics == 1) { 92 $scale *= 0.2; 93 $ntics = int($range/$scale); 94 } elsif ($ntics == 2) { 95 $scale *= 0.4; 96 $ntics = int($range/$scale); 97 } elsif ($ntics == 3) { 98 $scale *= 0.4; 99 $ntics = int($range/$scale); 100 } elsif ($ntics == 4) { 101 $scale *= 0.5; 102 $ntics = int($range/$scale); 103 } elsif ($ntics == 5) { 104 $scale *= 0.8; 105 $ntics = int($range/$scale); 106 } 107 my $lowval = ceil($low/$scale)*$scale; 108 my $highval = floor($high/$scale)*$scale; 109 my $ticstring = "$lowval,$scale,$highval"; 110 return $ticstring; 111} 112 113 114sub getMinMax { 115 my $ArrayRef = shift; 116 my $column = shift; 117 my $start = shift; 118 my $stop = shift; 119 120 if ($stop < 0) { 121 $stop = $#{$ArrayRef}; 122 } 123 124 my $min = 1000000000; 125 my $max = -1000000000; 126 my $numsamples = 0; 127 for (my $i = $start+1; $i <= $stop; $i++) { 128 my $data = ${$ArrayRef}[$i][$column]; 129 if ($data < $min) { 130 $min = $data; 131 } 132 if ($data > $max) { 133 $max = $data; 134 } 135 $numsamples++; 136 } 137 return ($min, $max, $numsamples); 138} 139 140sub getGPLHistString { 141 my $ArrayRef = shift; 142 my $column = shift; 143 my $start = shift; 144 my $stop = shift; 145 146 my $min; 147 my $max; 148 my $numsamples; 149 ($min, $max, $numsamples) = getMinMax($ArrayRef,$column,$start,$stop); 150 my @arr = getColumn($ArrayRef,$column,$start,$stop); 151 152 my $numbins = floor($numsamples / 50); 153 154 my @histdata; 155 for (my $i = 0; $i < $numbins; $i++) { 156 push @histdata, 0; 157 } 158 159 my $binwidth = ($max-$min)/($numbins+0.0); 160 my $maxcount = 0; 161 foreach my $i (@arr) { 162 my $binnum = floor(($i-$min)/$binwidth); 163 if ($binnum == $numbins || $binnum < 0) { 164 $binnum--; 165 } 166 $histdata[$binnum]++; 167 if ($histdata[$binnum] > $maxcount) { 168 $maxcount = $histdata[$binnum]; 169 } 170 } 171 my $plotmax = $maxcount+2; 172 173# my $plotstring = "unset ytics; set y2tics\n"; 174 my $plotstring .= "plot [$plotmax:0][$min:$max] \"-\" u 1:2 w l notitle\n "; 175# my $plotstring .= "plot [10:0][$min:$max] \"-\" u 1:2 w l notitle\n "; 176 for (my $i = 0; $i < $numbins; $i++) { 177 my $boxxmin = $min+$i*$binwidth; 178 my $boxxmax = $min+($i+1)*$binwidth; 179 my $boxymax = $histdata[$i]; 180 my $tag = $i+1; 181 $plotstring .= "0 $boxxmin\n $boxymax $boxxmin\n $boxymax $boxxmax\n 0 $boxxmax\n"; 182 } 183 $plotstring .= "end \n"; 184 return $plotstring; 185} 186 187sub getColumn { 188 my $matrix = shift; 189 my $colnum = shift; 190 my $startnum = shift; 191 my $end = shift; 192 193 if ($end < 0) { 194 $end = $#{$matrix}; 195 } 196 197 $startnum || ($startnum = 0); 198 my @outmat; 199 for (my $i = $startnum+1; $i <= $end; $i++) { 200 push @outmat, ${$matrix}[$i][$colnum]; 201 } 202 return @outmat; 203} 204 205 206 207