1#! /usr/bin/env python3
2from __future__ import print_function
3
4# Statical error checking code for use by testing framework
5# Jaron Krogel/ORNL
6
7# To maximize portability, only standard Python modules should
8# be used (that is, no numpy)
9
10import os
11import sys
12from optparse import OptionParser
13import math
14
15
16# standalone definition of error function from Abramowitz & Stegun
17# credit: http://www.johndcook.com/blog/2009/01/19/stand-alone-error-function-erf/
18def erf(x):
19    # constants
20    a1 =  0.254829592
21    a2 = -0.284496736
22    a3 =  1.421413741
23    a4 = -1.453152027
24    a5 =  1.061405429
25    p  =  0.3275911
26
27    # Save the sign of x
28    sign = 1
29    if x < 0:
30        sign = -1
31    x = abs(x)
32
33    # A & S 7.1.26
34    t = 1.0/(1.0 + p*x)
35    y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*math.exp(-x*x)
36
37    return sign*y
38#end def erf
39
40
41# Returns failure error code to OS.
42# Explicitly prints 'fail' after an optional message.
43def exit_fail(msg=None):
44    if msg!=None:
45        print(msg)
46    #end if
47    print('Test status: fail')
48    exit(1)
49#end def exit_fail
50
51
52# Returns success error code to OS.
53# Explicitly prints 'pass' after an optional message.
54def exit_pass(msg=None):
55    if msg!=None:
56        print(msg)
57    #end if
58    print('Test status: pass')
59    exit(0)
60#end def exit_pass
61
62def compute_mean(v):
63    if len(v) == 0:
64        return 0.0
65    return sum(v)/len(v)
66
67
68def compute_variance(v):
69    if len(v) == 0:
70        return 0.0
71    mean = compute_mean(v)
72    return sum([(x-mean)**2 for x in v])/len(v)
73
74
75
76# Calculates the mean, variance, errorbar, and autocorrelation time
77# for a 1-d array of statistical data values.
78# If 'exclude' is provided, the first 'exclude' values will be
79# excluded from the analysis.
80def simstats(x,exclude=None):
81    if exclude!=None:
82        x = x[exclude:]
83    #end if
84
85    N    = len(x)
86    mean = compute_mean(x)
87    var = compute_variance(x)
88
89    i=0
90    tempC=0.5
91    kappa=0.0
92    if abs(var)<1e-15:
93        kappa = 1.0
94    else:
95        ovar=1.0/var
96        while (tempC>0 and i<(N-1)):
97            kappa=kappa+2.0*tempC
98            i=i+1
99            tempC = ovar/(N-i)*sum([ (x[idx]-mean)*(x[idx+i]-mean) for idx in range(N-i) ])
100        #end while
101        if kappa == 0.0:
102            kappa = 1.0
103        #end if
104    #end if
105    Neff=(N+0.0)/(kappa+0.0)
106    if (Neff == 0.0):
107        Neff = 1.0
108    #end if
109    error=math.sqrt(var/Neff)
110
111    return (mean,var,error,kappa)
112#end def simstats
113
114
115
116# Reads command line options.
117# For example:
118#   check_scalars.py  --ns 2  -p Li  -s '1 3' -e 10  --le '-7.478011 0.000035  -7.478059 0.000035'
119# This invocation expects two scalar files to be present:
120#   Li.s001.scalar.dat  Li.s003.scalar.dat
121# The local energy ('le') will be computed excluding ('e') 10 blocks each.
122# the test passes if both of the following are true:
123#   |local_energy_mean_of_series_1 - (-7.478011)| < 2*0.000035
124#   |local_energy_mean_of_series_3 - (-7.478059)| < 2*0.000035
125# Note that the factor of 2 is specified by 'ns', the allowed number of sigma deviations for the test.
126# The inputted errorbar (sigma) to check against (0.000035) should satisfy the following:
127#   Let err_ref be the error bar of the reference solution.
128#   Let err_comp be the expected errorbar of the completed test run.
129#   Then the provided errorbar/sigma should be sqrt( err_ref**2 + err_comp**2 ).
130def read_command_line():
131    try:
132
133        parser = OptionParser(
134            usage='usage: %prog [options]',
135            add_help_option=False,
136            version='%prog 0.1'
137            )
138
139        parser.add_option('-h','--help',dest='help',
140                          action='store_true',default=False,
141                          help='Print help information and exit (default=%default).'
142                          )
143        parser.add_option('-p','--prefix',dest='prefix',
144                          default='qmc',
145                          help='Prefix for output files (default=%default).'
146                          )
147        parser.add_option('-s','--series',dest='series',
148                          default='0',
149                          help='Output series to analyze (default=%default).'
150                          )
151        parser.add_option('-e','--equilibration',dest='equilibration',
152                          default='0',
153                          help='Equilibration length in blocks (default=%default).'
154                          )
155        parser.add_option('-n','--nsigma',dest='nsigma',
156                          default='3',
157                          help='Sigma requirement for pass/fail (default=%default).'
158                          )
159
160        quantities = dict(
161            ar   = 'AcceptRatio',
162            le   = 'LocalEnergy',
163            va   = 'Variance',
164            ke   = 'Kinetic',
165            lp   = 'LocalPotential',
166            ee   = 'ElecElec',
167            cl   = 'Coulomb',
168            ii   = 'IonIon',
169            lpp  = 'LocalECP',
170            nlpp = 'NonLocalECP',
171            mpc  = 'MPC',
172            kec  = 'KEcorr',
173            bw   = 'BlockWeight',
174            ts   = 'TotalSamples',
175            fl   = 'Flux',
176            latdev = 'latdev',
177#now for some RMC estimators
178            ke_m = "Kinetic_m",
179            ke_p = "Kinetic_p",
180            ee_m = "ElecElec_m",
181            ee_p = "ElecElec_p",
182            lp_p = "LocalPotential_pure",
183#and some CSVMC estimators
184            le_A = "LocEne_0",
185            le_B = "LocEne_1",
186            dle_AB = "dLocEne_0_1",
187            ii_A = "IonIon_0",
188            ii_B = "IonIon_1",
189            dii_AB = "dIonIon_0_1",
190            ee_A = "ElecElec_0",
191            ee_B = "ElecElec_1",
192            dee_AB = "dElecElec_0_1",
193#AFQMC quantities
194            eloc    = 'Eloc',
195            elocest = 'ElocEstim',
196            el = 'EnergyEstim__nume_real',
197            )
198
199        for qshort in sorted(quantities.keys()):
200            qlong = quantities[qshort]
201            parser.add_option('--'+qshort,'--'+qlong,dest=qlong,
202                              default=None,
203                              help='Reference value and errorbar for '+qlong+' (one value/error pair per series).'
204                              )
205        #end for
206
207        # Check arbitrary named values in the scalar.dat file
208        parser.add_option('--name',default=None, help='Name of column to check in the scalar.dat file')
209        parser.add_option('--ref-value',default=None, help='Reference value for <name>')
210        parser.add_option('--ref-error',default=None, help='Reference error for <name>')
211
212        options,files_in = parser.parse_args()
213
214        if options.help:
215            print('\n'+parser.format_help().strip())
216            exit()
217        #end if
218
219        options.series         = [int(a) for a in options.series.split()]
220        options.equilibration  = [int(a) for a in options.equilibration.split()]
221        if len(options.series)>0 and len(options.equilibration)==1:
222            options.equilibration = len(options.series)*[options.equilibration[0]]
223        #end if
224        options.nsigma         = float(options.nsigma)
225
226        quants_check = []
227        for q in quantities.values():
228            v = options.__dict__[q]
229            if v!=None:
230                vref = [float(a) for a in v.split()]
231                if len(vref)!=2*len(options.series):
232                    exit_fail('must provide one reference value and errorbar for '+q)
233                #end if
234                options.__dict__[q] = vref
235                quants_check.append(q)
236            #end if
237        #end for
238
239
240        if options.name is not None and (options.ref_value is None or options.ref_error is None):
241            exit_fail("All of --name,--ref-value, and --ref-error options must be present")
242
243        if options.name and options.ref_value:
244            val = float(options.ref_value)
245            if options.ref_error:
246                err = float(options.ref_error)
247                options.__dict__[options.name] = [val,err]
248                quants_check.append(options.name)
249            #end if
250        #end if
251    except Exception as e:
252        exit_fail('error during command line read:\n'+str(e))
253    #end try
254
255    if len(quants_check)==0:
256        cmd = ''
257        for arg in sys.argv:
258            cmd += arg+' '
259        #end for
260        exit_fail('no quantities requested to check\nfull command: '+cmd)
261    #end if
262
263    return options,quants_check
264#end def read_command_line
265
266
267
268# Reads scalar.dat files and performs statistical analysis.
269def process_scalar_files(options,quants_check):
270    values = dict()
271
272    try:
273        ns = 0
274        for s in options.series:
275            svals = dict()
276            scalar_file = options.prefix+'.s'+str(int(s)).zfill(3)+'.scalar.dat'
277
278            if os.path.exists(scalar_file):
279                fobj = open(scalar_file,'r')
280                quantities = fobj.readline().split()[2:]
281                rawdata_list = []
282                for line in fobj:
283                    vals = line.strip().split()[1:]
284                    fp_vals = [float(v) for v in vals]
285                    rawdata_list.append(fp_vals)
286                fobj.close()
287
288                rawdata = [[]]
289                if len(rawdata_list) == 0:
290                    exit_fail('scalar file has no data: '+scalar_file)
291                # end if
292
293                # transpose
294                ncols = len(rawdata_list[0])
295                rawdata = [[] for i in range(ncols)]
296                for line in rawdata_list:
297                    for i in range(ncols):
298                        rawdata[i].append(line[i])
299                    # end for
300                # end for
301
302
303                equil = options.equilibration[ns]
304
305                data  = dict()
306                stats = dict()
307                for i in range(len(quantities)):
308                    q = quantities[i]
309                    d = rawdata[i]
310                    data[q]  = d
311                    stats[q] = simstats(d,equil)
312                #end for
313
314                if 'LocalEnergy_sq' in data and 'LocalEnergy' in data:
315                    v = [(le_sq-le**2) for le_sq,le in zip(data['LocalEnergy_sq'],data['LocalEnergy'])]
316                    stats['Variance'] = simstats(v,equil)
317                #end if
318
319                if 'BlockWeight' in data:
320                    ts = sum(data['BlockWeight'])
321                    stats['TotalSamples'] = (ts,0.0,0.0,1.0) # mean, var, error, kappa
322                #end if
323
324                for q in quants_check:
325                    if q in stats:
326                        mean,var,error,kappa = stats[q]
327                        svals[q] = mean,error
328                    else:
329                        exit_fail('{0} is not present in file {1}'.format(q,scalar_file))
330                    #end if
331                #end for
332            else:
333                exit_fail('scalar file does not exist: '+scalar_file)
334            #end if
335
336            values[s] = svals
337            ns += 1
338        #end for
339    except Exception as e:
340        exit_fail('error during scalar file processing:\n'+str(e))
341    #end try
342
343    return values
344#end def process_scalar_files
345
346
347# Checks computed values from scalar.dat files
348# against specified reference values.
349passfail = {True:'pass',False:'fail'}
350def check_values(options,quants_check,values):
351    success = True
352    msg = ''
353
354    try:
355        ns = 0
356        for s in options.series:
357            msg+='Tests for series {0}\n'.format(s)
358            for q in quants_check:
359                msg+='  Testing quantity: {0}\n'.format(q)
360
361                ref = options.__dict__[q]
362                mean_ref  = ref[2*ns]
363                error_ref = ref[2*ns+1]
364                mean_comp,error_comp = values[s][q]
365
366                quant_success = abs(mean_comp-mean_ref) <= options.nsigma*error_ref
367
368                success &= quant_success
369
370                delta = mean_comp-mean_ref
371                delta_err = math.sqrt(error_comp**2+error_ref**2)
372
373                msg+='    reference mean value     : {0: 12.8f}\n'.format(mean_ref)
374                msg+='    reference error bar      : {0: 12.8f}\n'.format(error_ref)
375                msg+='    computed  mean value     : {0: 12.8f}\n'.format(mean_comp)
376                msg+='    computed  error bar      : {0: 12.8f}\n'.format(error_comp)
377                msg+='    pass tolerance           : {0: 12.8f}  ({1: 12.8f} sigma)\n'.format(options.nsigma*error_ref,options.nsigma)
378                if error_ref > 0.0:
379                    msg+='    deviation from reference : {0: 12.8f}  ({1: 12.8f} sigma)\n'.format(delta,delta/error_ref)
380                msg+='    error bar of deviation   : {0: 12.8f}\n'.format(delta_err)
381                if error_ref > 0.0:
382                    msg+='    significance probability : {0: 12.8f}  (gaussian statistics)\n'.format(erf(abs(delta/error_ref)/math.sqrt(2.0)))
383                msg+='    status of this test      :   {0}\n'.format(passfail[quant_success])
384            #end for
385            ns+=1
386        #end for
387    except Exception as e:
388        exit_fail('error during value check:\n'+str(e))
389    #end try
390
391    return success,msg
392#end def check_values
393
394
395
396# Main execution
397if __name__=='__main__':
398    # Read and interpret command line options.
399    options,quants_check = read_command_line()
400
401    # Compute means of desired quantities from scalar.dat files.
402    values = process_scalar_files(options,quants_check)
403
404    # Check computed means against reference solutions.
405    success,msg = check_values(options,quants_check,values)
406
407    # Pass success/failure exit codes and strings to the OS.
408    if success:
409        exit_pass(msg)
410    else:
411        exit_fail(msg)
412    #end if
413#end if
414