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