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