1#! /usr/bin/env python3 2 3from __future__ import print_function, division 4 5# Statical error checking code for use by testing framework (stat.h5 files) 6# Jaron Krogel/ORNL 7 8 9# check_stats.py packages obj and HDFreader classes from Nexus. 10# Note that h5py is required (which depends on numpy). 11 12 13 14###################################################################### 15# from generic.py 16###################################################################### 17 18import sys 19import traceback 20from copy import deepcopy 21import pickle 22from random import randint 23 24 25class generic_settings: 26 devlog = sys.stdout 27 raise_error = False 28#end class generic_settings 29 30 31class NexusError(Exception): 32 None 33#end class NexusError 34 35 36exit_call = sys.exit 37 38 39def nocopy(value): 40 return value 41#end def nocopy 42 43 44 45def log(*items,**kwargs): 46 indent=None 47 logfile=generic_settings.devlog 48 if len(kwargs)>0: 49 n=0 50 if 'indent' in kwargs: 51 indent = kwargs['indent'] 52 n+=1 53 #end if 54 if 'logfile' in kwargs: 55 logfile = kwargs['logfile'] 56 n+=1 57 #end if 58 if n!=len(kwargs): 59 valid = 'indent logfile'.split() 60 invalid = set(kwargs.keys())-set(valid) 61 error('invalid keyword arguments provided\ninvalid keywords: {0}\nvalid options are: {1}'.format(sorted(invalid),valid)) 62 #end if 63 #end if 64 if len(items)==1 and isinstance(items[0],str): 65 s = items[0] 66 else: 67 s='' 68 for item in items: 69 s+=str(item)+' ' 70 #end for 71 #end if 72 if len(s)>0: 73 if isinstance(indent,str): 74 s=indent+s.replace('\n','\n'+indent) 75 #end if 76 s += '\n' 77 #end if 78 logfile.write(s) 79#end def log 80 81 82def message(msg,header=None,post_header=' message:',indent=' ',logfile=None): 83 if logfile is None: 84 logfile = generic_settings.devlog 85 #end if 86 if header is None: 87 header = post_header.lstrip() 88 else: 89 header += post_header 90 #end if 91 log('\n '+header,logfile=logfile) 92 log(msg.rstrip(),indent=indent,logfile=logfile) 93#end def message 94 95 96def warn(msg,header=None,indent=' ',logfile=None): 97 if logfile is None: 98 logfile = generic_settings.devlog 99 #end if 100 post_header=' warning:' 101 message(msg,header,post_header,indent,logfile) 102#end def warn 103 104 105def error(msg,header=None,exit=True,trace=True,indent=' ',logfile=None): 106 if generic_settings.raise_error: 107 raise NexusError(msg) 108 #end if 109 if logfile is None: 110 logfile = generic_settings.devlog 111 #end if 112 post_header=' error:' 113 message(msg,header,post_header,indent,logfile) 114 if exit: 115 log(' exiting.\n') 116 if trace: 117 traceback.print_stack() 118 #end if 119 exit_call() 120 #end if 121#end def error 122 123 124 125class object_interface(object): 126 _logfile = sys.stdout 127 128 def __len__(self): 129 return len(self.__dict__) 130 #end def __len__ 131 132 def __contains__(self,name): 133 return name in self.__dict__ 134 #end def 135 136 def __getitem__(self,name): 137 return self.__dict__[name] 138 #end def __getitem__ 139 140 def __setitem__(self,name,value): 141 self.__dict__[name]=value 142 #end def __setitem__ 143 144 def __delitem__(self,name): 145 del self.__dict__[name] 146 #end def __delitem__ 147 148 def __iter__(self): 149 for item in self.__dict__: 150 yield self.__dict__[item] 151 #end for 152 #end def __iter__ 153 154 def __repr__(self): 155 s='' 156 for k in sorted(self._keys()): 157 if not isinstance(k,str) or k[0]!='_': 158 v=self.__dict__[k] 159 if hasattr(v,'__class__'): 160 s+=' {0:<20} {1:<20}\n'.format(k,v.__class__.__name__) 161 else: 162 s+=' {0:<20} {1:<20}\n'.format(k,type(v)) 163 #end if 164 #end if 165 #end for 166 return s 167 #end def __repr__ 168 169 def __str__(self,nindent=1): 170 pad = ' ' 171 npad = nindent*pad 172 s='' 173 normal = [] 174 qable = [] 175 for k,v in self._items(): 176 if not isinstance(k,str) or k[0]!='_': 177 if isinstance(v,object_interface): 178 qable.append(k) 179 else: 180 normal.append(k) 181 #end if 182 #end if 183 #end for 184 normal.sort() 185 qable.sort() 186 indent = npad+18*' ' 187 for k in normal: 188 v = self[k] 189 vstr = str(v).replace('\n','\n'+indent) 190 s+=npad+'{0:<15} = '.format(k)+vstr+'\n' 191 #end for 192 for k in qable: 193 v = self[k] 194 s+=npad+str(k)+'\n' 195 s+=v.__str__(nindent+1) 196 if isinstance(k,str): 197 s+=npad+'end '+k+'\n' 198 #end if 199 #end for 200 return s 201 #end def __str__ 202 203 def __eq__(self,other): 204 if not hasattr(other,'__dict__'): 205 return False 206 #end if 207 eq = True 208 for sname in self.__dict__: 209 if sname not in other.__dict__: 210 return False 211 #end if 212 svar = self.__dict__[sname] 213 ovar = other.__dict__[sname] 214 stype = type(svar) 215 otype = type(ovar) 216 if stype!=otype: 217 return False 218 #end if 219 eqval = svar==ovar 220 if isinstance(eqval,bool): 221 eq &= eqval 222 else: 223 try: # accommodate numpy arrays implicitly 224 eq &= eqval.all() 225 except: 226 return False 227 #end try 228 #end if 229 #end for 230 return eq 231 #end def __eq__ 232 233 def tree(self,depth=None,all=False,types=False,nindent=1): 234 if depth==nindent-1: 235 return '' 236 #end if 237 pad = ' ' 238 npad = nindent*pad 239 s='' 240 normal = [] 241 qable = [] 242 for k,v in self._items(): 243 if not isinstance(k,str) or k[0]!='_': 244 if isinstance(v,object_interface): 245 qable.append(k) 246 else: 247 normal.append(k) 248 #end if 249 #end if 250 #end for 251 normal.sort() 252 qable.sort() 253 indent = npad+18*' ' 254 if all: 255 for k in normal: 256 v = self[k] 257 if types: 258 s+=npad+'{0:<15} = '.format(k) 259 if hasattr(v,'__class__'): 260 s+='{0:<20}'.format(v.__class__.__name__) 261 else: 262 s+='{0:<20}'.format(type(v)) 263 #end if 264 else: 265 s+=npad+str(k) 266 #end if 267 s+='\n' 268 #end for 269 #end if 270 if all and depth!=nindent: 271 for k in qable: 272 v = self[k] 273 s+=npad+str(k)+'\n' 274 s+=v.tree(depth,all,types,nindent+1) 275 if isinstance(k,str): 276 s+=npad+'end '+k+'\n' 277 #end if 278 #end for 279 else: 280 for k in qable: 281 v = self[k] 282 if types: 283 s+=npad+'{0:<15} = '.format(k) 284 if hasattr(v,'__class__'): 285 s+='{0:<20}'.format(v.__class__.__name__) 286 else: 287 s+='{0:<20}'.format(type(v)) 288 #end if 289 else: 290 s+=npad+str(k) 291 #end if 292 s+='\n' 293 s+=v.tree(depth,all,types,nindent+1) 294 #end for 295 #end if 296 return s 297 #end def tree 298 299 300 # dict interface 301 def keys(self): 302 return self.__dict__.keys() 303 #end def keys 304 305 def values(self): 306 return self.__dict__.values() 307 #end def values 308 309 def items(self): 310 return self.__dict__.items() 311 #end def items 312 313 def copy(self): 314 return deepcopy(self) 315 #end def copy 316 317 def clear(self): 318 self.__dict__.clear() 319 #end def clear 320 321 322 # save/load 323 def save(self,fpath=None): 324 if fpath==None: 325 fpath='./'+self.__class__.__name__+'.p' 326 #end if 327 fobj = open(fpath,'w') 328 binary = pickle.HIGHEST_PROTOCOL 329 pickle.dump(self,fobj,binary) 330 fobj.close() 331 del fobj 332 del binary 333 return 334 #end def save 335 336 def load(self,fpath=None): 337 if fpath==None: 338 fpath='./'+self.__class__.__name__+'.p' 339 #end if 340 fobj = open(fpath,'r') 341 tmp = pickle.load(fobj) 342 fobj.close() 343 d = self.__dict__ 344 d.clear() 345 for k,v in tmp.__dict__.items(): 346 d[k] = v 347 #end for 348 del fobj 349 del tmp 350 return 351 #end def load 352 353 354 355 # log, warning, and error messages 356 def open_log(self,filepath): 357 self._logfile = open(filepath,'w') 358 #end def open_log 359 360 def close_log(self): 361 self._logfile.close() 362 #end def close_log 363 364 def write(self,s): 365 self._logfile.write(s) 366 #end def write 367 368 def log(self,*items,**kwargs): 369 if 'logfile' not in kwargs: 370 kwargs['logfile'] = self._logfile 371 #end if 372 log(*items,**kwargs) 373 #end def log 374 375 def warn(self,message,header=None): 376 if header is None: 377 header=self.__class__.__name__ 378 #end if 379 warn(message,header,logfile=self._logfile) 380 #end def warn 381 382 def error(self,message,header=None,exit=True,trace=True): 383 if header==None: 384 header = self.__class__.__name__ 385 #end if 386 error(message,header,exit,trace,logfile=self._logfile) 387 #end def error 388 389 @classmethod 390 def class_log(cls,message): 391 log(message,logfile=cls._logfile) 392 #end def class_log 393 394 @classmethod 395 def class_warn(cls,message,header=None,post_header=' Warning:'): 396 if header==None: 397 header=cls.__name__ 398 #end if 399 warn(message,header,logfile=cls._logfile) 400 #end def class_warn 401 402 @classmethod 403 def class_error(cls,message,header=None,exit=True,trace=True,post_header=' Error:'): 404 if header==None: 405 header = cls.__name__ 406 #end if 407 error(message,header,exit,trace,logfile=cls._logfile) 408 #end def class_error 409 410 @classmethod 411 def class_has(cls,k): 412 return hasattr(cls,k) 413 #end def classmethod 414 415 @classmethod 416 def class_keys(cls): 417 return cls.__dict__.keys() 418 #end def class_keys 419 420 @classmethod 421 def class_get(cls,k): 422 return getattr(cls,k) 423 #end def class_set 424 425 @classmethod 426 def class_set(cls,**kwargs): 427 for k,v in kwargs.items(): 428 setattr(cls,k,v) 429 #end for 430 #end def class_set 431 432 @classmethod 433 def class_set_single(cls,k,v): 434 setattr(cls,k,v) 435 #end def class_set_single 436 437 @classmethod 438 def class_set_optional(cls,**kwargs): 439 for k,v in kwargs.items(): 440 if not hasattr(cls,k): 441 setattr(cls,k,v) 442 #end if 443 #end for 444 #end def class_set_optional 445 446 447 # access preserving functions 448 # dict interface 449 def _keys(self,*args,**kwargs): 450 return object_interface.keys(self,*args,**kwargs) 451 def _values(self,*args,**kwargs): 452 object_interface.values(self,*args,**kwargs) 453 def _items(self,*args,**kwargs): 454 return object_interface.items(self,*args,**kwargs) 455 def _copy(self,*args,**kwargs): 456 return object_interface.copy(self,*args,**kwargs) 457 def _clear(self,*args,**kwargs): 458 object_interface.clear(self,*args,**kwargs) 459 # save/load 460 def _save(self,*args,**kwargs): 461 object_interface.save(self,*args,**kwargs) 462 def _load(self,*args,**kwargs): 463 object_interface.load(self,*args,**kwargs) 464 # log, warning, and error messages 465 def _open_log(self,*args,**kwargs): 466 object_interface.open_log(self,*args,**kwargs) 467 def _close_log(self,*args,**kwargs): 468 object_interface.close_log(self,*args,**kwargs) 469 def _write(self,*args,**kwargs): 470 object_interface.write(self,*args,**kwargs) 471 def _log(self,*args,**kwargs): 472 object_interface.log(self,*args,**kwargs) 473 def _error(self,*args,**kwargs): 474 object_interface.error(self,*args,**kwargs) 475 def _warn(self,*args,**kwargs): 476 object_interface.warn(self,*args,**kwargs) 477 478#end class object_interface 479 480 481 482class obj(object_interface): 483 484 def __init__(self,*vars,**kwargs): 485 for var in vars: 486 if isinstance(var,(dict,object_interface)): 487 for k,v in var.items(): 488 self[k] = v 489 #end for 490 else: 491 self[var] = None 492 #end if 493 #end for 494 for k,v in kwargs.items(): 495 self[k] = v 496 #end for 497 #end def __init__ 498 499 500 # list interface 501 def append(self,value): 502 self[len(self)] = value 503 #end def append 504 505 506 # return representations 507 def list(self,*keys): 508 nkeys = len(keys) 509 if nkeys==0: 510 keys = sorted(self._keys()) 511 elif nkeys==1 and isinstance(keys[0],(list,tuple)): 512 keys = keys[0] 513 #end if 514 values = [] 515 for key in keys: 516 values.append(self[key]) 517 #end if 518 return values 519 #end def list 520 521 def list_optional(self,*keys): 522 nkeys = len(keys) 523 if nkeys==0: 524 keys = sorted(self._keys()) 525 elif nkeys==1 and isinstance(keys[0],(list,tuple)): 526 keys = keys[0] 527 #end if 528 values = [] 529 for key in keys: 530 if key in self: 531 values.append(self[key]) 532 else: 533 values.append(None) 534 #end if 535 #end if 536 return values 537 #end def list_optional 538 539 def tuple(self,*keys): 540 return tuple(obj.list(self,*keys)) 541 #end def tuple 542 543 def dict(self,*keys): 544 nkeys = len(keys) 545 if nkeys==0: 546 keys = sorted(self._keys()) 547 elif nkeys==1 and isinstance(keys[0],(list,tuple)): 548 keys = keys[0] 549 #end if 550 d = dict() 551 for k in keys: 552 d[k] = self[k] 553 #end for 554 return d 555 #end def dict 556 557 def to_dict(self): 558 d = dict() 559 for k,v in self._items(): 560 if isinstance(v,obj): 561 d[k] = v._to_dict() 562 else: 563 d[k] = v 564 #end if 565 #end for 566 return d 567 #end def to_dict 568 569 def obj(self,*keys): 570 nkeys = len(keys) 571 if nkeys==0: 572 keys = sorted(self._keys()) 573 elif nkeys==1 and isinstance(keys[0],(list,tuple)): 574 keys = keys[0] 575 #end if 576 o = obj() 577 for k in keys: 578 o[k] = self[k] 579 #end for 580 return o 581 #end def obj 582 583 584 # list extensions 585 def first(self): 586 return self[min(self._keys())] 587 #end def first 588 589 def last(self): 590 return self[max(self._keys())] 591 #end def last 592 593 def select_random(self): 594 return self[randint(0,len(self)-1)] 595 #end def select_random 596 597 598 # dict extensions 599 def random_key(self): 600 key = None 601 nkeys = len(self) 602 if nkeys>0: 603 key = self._keys()[randint(0,nkeys-1)] 604 #end if 605 return key 606 #end def random_key 607 608 609 def set(self,*objs,**kwargs): 610 for key,value in kwargs.items(): 611 self[key]=value 612 #end for 613 if len(objs)>0: 614 for o in objs: 615 for k,v in o.items(): 616 self[k] = v 617 #end for 618 #end for 619 #end if 620 return self 621 #end def set 622 623 def set_optional(self,*objs,**kwargs): 624 for key,value in kwargs.items(): 625 if key not in self: 626 self[key]=value 627 #end if 628 #end for 629 if len(objs)>0: 630 for o in objs: 631 for k,v in o.items(): 632 if k not in self: 633 self[k] = v 634 #end if 635 #end for 636 #end for 637 #end if 638 return self 639 #end def set_optional 640 641 def get(self,key,value=None): # follow dict interface, no plural 642 if key in self: 643 value = self[key] 644 #end if 645 return value 646 #end def get 647 648 def get_optional(self,key,value=None): 649 if key in self: 650 value = self[key] 651 #end if 652 return value 653 #end def get_optional 654 655 def get_required(self,key): 656 if key in self: 657 value = self[key] 658 else: 659 obj.error(self,'a required key is not present\nkey required: {0}\nkeys present: {1}'.format(key,sorted(self._keys()))) 660 #end if 661 return value 662 #end def get_required 663 664 def delete(self,*keys): 665 nkeys = len(keys) 666 single = False 667 if nkeys==0: 668 keys = sorted(self._keys()) 669 elif nkeys==1 and isinstance(keys[0],(list,tuple)): 670 keys = keys[0] 671 elif nkeys==1: 672 single = True 673 #end if 674 values = [] 675 for key in keys: 676 values.append(self[key]) 677 del self[key] 678 #end for 679 if single: 680 return values[0] 681 else: 682 return values 683 #end if 684 #end def delete 685 686 def delete_optional(self,key,value=None): 687 if key in self: 688 value = self[key] 689 del self[key] 690 #end if 691 return value 692 #end def delete_optional 693 694 def delete_required(self,key): 695 if key in self: 696 value = self[key] 697 del self[key] 698 else: 699 obj.error(self,'a required key is not present\nkey required: {0}\nkeys present: {1}'.format(key,sorted(self._keys()))) 700 #end if 701 return value 702 #end def delete_required 703 704 def add(self,key,value): 705 self[key] = value 706 #end def add 707 708 def add_optional(self,key,value): 709 if key not in self: 710 self[key] = value 711 #end if 712 #end def add_optional 713 714 def transfer_from(self,other,keys=None,copy=False,overwrite=True): 715 if keys==None: 716 if isinstance(other,object_interface): 717 keys = other._keys() 718 else: 719 keys = other.keys() 720 #end if 721 #end if 722 if copy: 723 copier = deepcopy 724 else: 725 copier = nocopy 726 #end if 727 if overwrite: 728 for k in keys: 729 self[k]=copier(other[k]) 730 #end for 731 else: 732 for k in keys: 733 if k not in self: 734 self[k]=copier(other[k]) 735 #end if 736 #end for 737 #end if 738 #end def transfer_from 739 740 def transfer_to(self,other,keys=None,copy=False,overwrite=True): 741 if keys==None: 742 keys = self._keys() 743 #end if 744 if copy: 745 copier = deepcopy 746 else: 747 copier = nocopy 748 #end if 749 if overwrite: 750 for k in keys: 751 other[k]=copier(self[k]) 752 #end for 753 else: 754 for k in keys: 755 if k not in self: 756 other[k]=copier(self[k]) 757 #end if 758 #end for 759 #end if 760 #end def transfer_to 761 762 def move_from(self,other,keys=None): 763 if keys==None: 764 if isinstance(other,object_interface): 765 keys = other._keys() 766 else: 767 keys = other.keys() 768 #end if 769 #end if 770 for k in keys: 771 self[k]=other[k] 772 del other[k] 773 #end for 774 #end def move_from 775 776 def move_to(self,other,keys=None): 777 if keys==None: 778 keys = self._keys() 779 #end if 780 for k in keys: 781 other[k]=self[k] 782 del self[k] 783 #end for 784 #end def move_to 785 786 def copy_from(self,other,keys=None,deep=True): 787 obj.transfer_from(self,other,keys,copy=deep) 788 #end def copy_from 789 790 def copy_to(self,other,keys=None,deep=True): 791 obj.transfer_to(self,other,keys,copy=deep) 792 #end def copy_to 793 794 def shallow_copy(self): 795 new = self.__class__() 796 for k,v in self._items(): 797 new[k] = v 798 #end for 799 return new 800 #end def shallow_copy 801 802 def inverse(self): 803 new = self.__class__() 804 for k,v in self._items(): 805 new[v] = k 806 #end for 807 return new 808 #end def inverse 809 810 def path_exists(self,path): 811 o = self 812 if isinstance(path,str): 813 path = path.split('/') 814 #end if 815 for p in path: 816 if not p in o: 817 return False 818 #end if 819 o = o[p] 820 #end for 821 return True 822 #end def path_exists 823 824 def set_path(self,path,value=None): 825 o = self 826 cls = self.__class__ 827 if isinstance(path,str): 828 path = path.split('/') 829 #end if 830 for p in path[0:-1]: 831 if not p in o: 832 o[p] = cls() 833 #end if 834 o = o[p] 835 #end for 836 o[path[-1]] = value 837 #end def set_path 838 839 def get_path(self,path,value=None): 840 o = self 841 if isinstance(path,str): 842 path = path.split('/') 843 #end if 844 for p in path[0:-1]: 845 if not p in o: 846 return value 847 #end if 848 o = o[p] 849 #end for 850 lp = path[-1] 851 if lp not in o: 852 return value 853 else: 854 return o[lp] 855 #end if 856 #end def get_path 857 858 def serial(self,s=None,path=None): 859 first = s is None 860 if first: 861 s = obj() 862 path = '' 863 #end if 864 for k,v in self._items(): 865 p = path+str(k) 866 if isinstance(v,obj): 867 if len(v)==0: 868 s[p]=v 869 else: 870 v._serial(s,p+'/') 871 #end if 872 else: 873 s[p]=v 874 #end if 875 #end for 876 if first: 877 return s 878 #end if 879 #end def serial 880 881 882 883 # access preserving functions 884 # list interface 885 def _append(self,*args,**kwargs): 886 obj.append(self,*args,**kwargs) 887 # return representations 888 def _list(self,*args,**kwargs): 889 return obj.list(self,*args,**kwargs) 890 def _list_optional(self,*args,**kwargs): 891 return obj.list_optional(self,*args,**kwargs) 892 def _tuple(self,*args,**kwargs): 893 return obj.tuple(self,*args,**kwargs) 894 def _dict(self,*args,**kwargs): 895 return obj.dict(self,*args,**kwargs) 896 def _to_dict(self,*args,**kwargs): 897 return obj.to_dict(self,*args,**kwargs) 898 def _obj(self,*args,**kwargs): 899 return obj.obj(self,*args,**kwargs) 900 # list extensions 901 def _first(self,*args,**kwargs): 902 return obj.first(self,*args,**kwargs) 903 def _last(self,*args,**kwargs): 904 return obj.last(self,*args,**kwargs) 905 def _select_random(self,*args,**kwargs): 906 return obj.select_random(self,*args,**kwargs) 907 # dict extensions 908 def _random_key(self,*args,**kwargs): 909 obj.random_key(self,*args,**kwargs) 910 def _set(self,*args,**kwargs): 911 obj.set(self,*args,**kwargs) 912 def _set_optional(self,*args,**kwargs): 913 obj.set_optional(self,*args,**kwargs) 914 def _get(self,*args,**kwargs): 915 obj.get(self,*args,**kwargs) 916 def _get_optional(self,*args,**kwargs): 917 obj.get_optional(self,*args,**kwargs) 918 def _get_required(self,*args,**kwargs): 919 obj.get_required(self,*args,**kwargs) 920 def _delete(self,*args,**kwargs): 921 obj.delete(self,*args,**kwargs) 922 def _delete_optional(self,*args,**kwargs): 923 obj.delete_optional(self,*args,**kwargs) 924 def _delete_required(self,*args,**kwargs): 925 obj.delete_required(self,*args,**kwargs) 926 def _add(self,*args,**kwargs): 927 obj.add(self,*args,**kwargs) 928 def _add_optional(self,*args,**kwargs): 929 obj.add_optional(self,*args,**kwargs) 930 def _transfer_from(self,*args,**kwargs): 931 obj.transfer_from(self,*args,**kwargs) 932 def _transfer_to(self,*args,**kwargs): 933 obj.transfer_to(self,*args,**kwargs) 934 def _move_from(self,*args,**kwargs): 935 obj.move_from(self,*args,**kwargs) 936 def _move_to(self,*args,**kwargs): 937 obj.move_to(self,*args,**kwargs) 938 def _copy_from(self,*args,**kwargs): 939 obj.copy_from(self,*args,**kwargs) 940 def _copy_to(self,*args,**kwargs): 941 obj.copy_to(self,*args,**kwargs) 942 def _shallow_copy(self,*args,**kwargs): 943 obj.shallow_copy(self,*args,**kwargs) 944 def _inverse(self,*args,**kwargs): 945 return obj.inverse(self,*args,**kwargs) 946 def _path_exists(self,*args,**kwargs): 947 obj.path_exists(self,*args,**kwargs) 948 def _set_path(self,*args,**kwargs): 949 obj.set_path(self,*args,**kwargs) 950 def _get_path(self,*args,**kwargs): 951 obj.get_path(self,*args,**kwargs) 952 def _serial(self,*args,**kwargs): 953 return obj.serial(self,*args,**kwargs) 954 955#end class obj 956 957###################################################################### 958# end from generic.py 959###################################################################### 960 961 962###################################################################### 963# from superstring.py 964###################################################################### 965 966import string 967 968def contains_any(str, set): 969 for c in set: 970 if c in str: return 1; 971 return 0; 972#end def contains_any 973 974invalid_variable_name_chars=set('!"#$%&\'()*+,-./:;<=>?@[\\]^`{|}-\n\t ') 975def valid_variable_name(s): 976 return not contains_any(s,invalid_variable_name_chars) 977#end def valid_variable_name 978 979###################################################################### 980# end from superstring.py 981###################################################################### 982 983 984###################################################################### 985# from debug.py 986###################################################################### 987 988import code 989import inspect 990 991def ci(locs=None,globs=None): 992 if locs is None or globs is None: 993 cur_frame = inspect.currentframe() 994 caller_frame = cur_frame.f_back 995 locs = caller_frame.f_locals 996 globs = caller_frame.f_globals 997 #end if 998 code.interact(local=dict(globs,**locs)) 999#end def ci 1000 1001ls = locals 1002gs = globals 1003interact = ci 1004 1005###################################################################### 1006# end from debug.py 1007###################################################################### 1008 1009 1010###################################################################### 1011# from developer.py 1012###################################################################### 1013 1014class DevBase(obj): 1015 def not_implemented(self): 1016 self.error('a base class function has not been implemented',trace=True) 1017 #end def not_implemented 1018#end class DevBase 1019 1020###################################################################### 1021# end from developer.py 1022###################################################################### 1023 1024 1025###################################################################### 1026# from hdfreader.py 1027###################################################################### 1028from numpy import array,ndarray,minimum,abs,ix_,resize 1029import sys 1030import keyword 1031from inspect import getmembers 1032import h5py 1033 1034 1035class HDFglobals(DevBase): 1036 view = False 1037#end class HDFglobals 1038 1039 1040class HDFgroup(DevBase): 1041 def _escape_name(self,name): 1042 if name in self._escape_names: 1043 name=name+'_' 1044 #end if 1045 return name 1046 #end def escape_name 1047 1048 def _set_parent(self,parent): 1049 self._parent=parent 1050 return 1051 #end def set_parent 1052 1053 def _add_dataset(self,name,dataset): 1054 self._datasets[name]=dataset 1055 return 1056 #end def add_dataset 1057 1058 def _add_group(self,name,group): 1059 group._name=name 1060 self._groups[name]=group 1061 return 1062 #end def add_group 1063 1064 def _contains_group(self,name): 1065 return name in self._groups.keys() 1066 #end def _contains_group 1067 1068 def _contains_dataset(self,name): 1069 return name in self._datasets.keys() 1070 #end def _contains_dataset 1071 1072 def _to_string(self): 1073 s='' 1074 if len(self._datasets)>0: 1075 s+=' datasets:\n' 1076 for k,v in self._datasets.items(): 1077 s+= ' '+k+'\n' 1078 #end for 1079 #end if 1080 if len(self._groups)>0: 1081 s+= ' groups:\n' 1082 for k,v in self._groups.items(): 1083 s+= ' '+k+'\n' 1084 #end for 1085 #end if 1086 return s 1087 #end def list 1088 1089# def __str__(self): 1090# return self._to_string() 1091# #end def __str__ 1092# 1093# def __repr__(self): 1094# return self._to_string() 1095# #end def __repr__ 1096 1097 def __init__(self): 1098 self._name='' 1099 self._parent=None 1100 self._groups={}; 1101 self._datasets={}; 1102 self._group_counts={} 1103 1104 self._escape_names=None 1105 self._escape_names=set(dict(getmembers(self)).keys()) | set(keyword.kwlist) 1106 return 1107 #end def __init__ 1108 1109 1110 def _remove_hidden(self,deep=True): 1111 if '_parent' in self: 1112 del self._parent 1113 #end if 1114 if deep: 1115 for name,value in self.items(): 1116 if isinstance(value,HDFgroup): 1117 value._remove_hidden() 1118 #end if 1119 #end for 1120 #end if 1121 for name in list(self.keys()): 1122 if name[0]=='_': 1123 del self[name] 1124 #end if 1125 #end for 1126 #end def _remove_hidden 1127 1128 1129 # read in all data views (h5py datasets) into arrays 1130 # useful for converting a single group read in view form to full arrays 1131 def read_arrays(self): 1132 self._remove_hidden() 1133 for k,v in self.items(): 1134 if isinstance(v,HDFgroup): 1135 v.read_arrays() 1136 else: 1137 self[k] = array(v) 1138 #end if 1139 #end for 1140 #end def read_arrays 1141 1142 1143 def get_keys(self): 1144 if '_groups' in self: 1145 keys = list(self._groups.keys()) 1146 else: 1147 keys = list(self.keys()) 1148 #end if 1149 return keys 1150 #end def get_keys 1151#end class HDFgroup 1152 1153 1154 1155 1156class HDFreader(DevBase): 1157 datasets = set(["<class 'h5py.highlevel.Dataset'>","<class 'h5py._hl.dataset.Dataset'>"]) 1158 groups = set(["<class 'h5py.highlevel.Group'>","<class 'h5py._hl.group.Group'>"]) 1159 1160 def __init__(self,fpath,verbose=False,view=False): 1161 1162 HDFglobals.view = view 1163 1164 if verbose: 1165 print(' Initializing HDFreader') 1166 1167 self.fpath=fpath 1168 if verbose: 1169 print(' loading h5 file') 1170 1171 try: 1172 self.hdf = h5py.File(fpath,'r') 1173 except IOError: 1174 self._success = False 1175 self.hdf = obj(obj=obj()) 1176 else: 1177 self._success = True 1178 #end if 1179 1180 if verbose: 1181 print(' converting h5 file to dynamic object') 1182 #convert the hdf 'dict' into a dynamic object 1183 self.nlevels=1 1184 self.ilevel=0 1185 # Set the current hdf group 1186 self.obj = HDFgroup() 1187 self.cur=[self.obj] 1188 self.hcur=[self.hdf] 1189 1190 if self._success: 1191 cur = self.cur[self.ilevel] 1192 hcur = self.hcur[self.ilevel] 1193 for kr,v in hcur.items(): 1194 k=cur._escape_name(kr) 1195 if valid_variable_name(k): 1196 vtype = str(type(v)) 1197 if vtype in HDFreader.datasets: 1198 self.add_dataset(cur,k,v) 1199 elif vtype in HDFreader.groups: 1200 self.add_group(hcur,cur,k,v) 1201 else: 1202 print('hdfreader error: encountered invalid type: '+vtype) 1203 sys.exit() 1204 #end if 1205 else: 1206 print('hdfreader warning: attribute '+k+' is not a valid variable name and has been ignored') 1207 #end if 1208 #end for 1209 #end if 1210 1211 if verbose: 1212 print(' end HDFreader Initialization') 1213 1214 return 1215 #end def __init__ 1216 1217 1218 def increment_level(self): 1219 self.ilevel+=1 1220 self.nlevels = max(self.ilevel+1,self.nlevels) 1221 if self.ilevel+1==self.nlevels: 1222 self.cur.append(None) 1223 self.hcur.append(None) 1224 #end if 1225 self.pad = self.ilevel*' ' 1226 return 1227 #end def increment_level 1228 1229 def decrement_level(self): 1230 self.ilevel-=1 1231 self.pad = self.ilevel*' ' 1232 return 1233 #end def decrement_level 1234 1235 def add_dataset(self,cur,k,v): 1236 if not HDFglobals.view: 1237 cur[k]=array(v) 1238 else: 1239 cur[k] = v 1240 #end if 1241 cur._add_dataset(k,cur[k]) 1242 return 1243 #end def add_dataset 1244 1245 def add_group(self,hcur,cur,k,v): 1246 cur[k] = HDFgroup() 1247 cur._add_group(k,cur[k]) 1248 cur._groups[k]._parent = cur 1249 self.increment_level() 1250 self.cur[self.ilevel] = cur._groups[k] 1251 self.hcur[self.ilevel] = hcur[k] 1252 1253 cur = self.cur[self.ilevel] 1254 hcur = self.hcur[self.ilevel] 1255 for kr,v in hcur.items(): 1256 k=cur._escape_name(kr) 1257 if valid_variable_name(k): 1258 vtype = str(type(v)) 1259 if vtype in HDFreader.datasets: 1260 self.add_dataset(cur,k,v) 1261 elif vtype in HDFreader.groups: 1262 self.add_group(hcur,cur,k,v) 1263 #end if 1264 else: 1265 print('hdfreader warning: attribute '+k+' is not a valid variable name and has been ignored') 1266 #end if 1267 #end for 1268 1269 return 1270 #end def add_group 1271#end class HDFreader 1272 1273 1274 1275def read_hdf(fpath,verbose=False,view=False): 1276 return HDFreader(fpath=fpath,verbose=verbose,view=view).obj 1277#end def read_hdf 1278 1279###################################################################### 1280# end from hdfreader.py 1281###################################################################### 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291import os 1292import sys 1293from optparse import OptionParser 1294from numpy import zeros,sqrt,longdouble,loadtxt 1295can_plot = False 1296try: 1297 import matplotlib 1298 gui_envs = ['GTKAgg','TKAgg','Qt4Agg','WXAgg'] 1299 for gui in gui_envs: 1300 try: 1301 matplotlib.use(gui,warn=False, force=True) 1302 from matplotlib import pyplot 1303 can_plot = True 1304 break 1305 except: 1306 continue 1307 #end try 1308 #end for 1309 from matplotlib.pyplot import figure,plot,xlabel,ylabel,title,show,ylim,legend,xlim,rcParams,savefig,bar,xticks,subplot,grid,setp,errorbar,loglog,semilogx,semilogy,text 1310 1311 params = {'legend.fontsize':14,'figure.facecolor':'white','figure.subplot.hspace':0., 1312 'axes.labelsize':16,'xtick.labelsize':14,'ytick.labelsize':14} 1313 rcParams.update(params) 1314except (ImportError,RuntimeError): 1315 can_plot = False 1316#end try 1317 1318 1319 1320class ColorWheel(DevBase): 1321 def __init__(self): 1322 colors = 'Black Maroon DarkOrange Green DarkBlue Purple Gray Firebrick Orange MediumSeaGreen DodgerBlue MediumOrchid'.split() 1323 lines = '- -- -. :'.split() 1324 markers = '. v s o ^ d p'.split() 1325 ls = [] 1326 for line in lines: 1327 for color in colors: 1328 ls.append((color,line)) 1329 #end for 1330 #end for 1331 ms = [] 1332 for i in range(len(markers)): 1333 ms.append((colors[i],markers[i])) 1334 #end for 1335 mls = [] 1336 ic=-1 1337 for line in lines: 1338 for marker in markers: 1339 ic = (ic+1)%len(colors) 1340 mls.append((colors[ic],marker+line)) 1341 #end for 1342 #end for 1343 self.line_styles = ls 1344 self.marker_styles = ms 1345 self.marker_line_styles = mls 1346 self.reset() 1347 #end def __init__ 1348 1349 def next_line(self): 1350 self.iline = (self.iline+1)%len(self.line_styles) 1351 return self.line_styles[self.iline] 1352 #end def next_line 1353 1354 def next_marker(self): 1355 self.imarker = (self.imarker+1)%len(self.marker_styles) 1356 return self.marker_styles[self.imarker] 1357 #end def next_marker 1358 1359 def next_marker_line(self): 1360 self.imarker_line = (self.imarker_line+1)%len(self.marker_line_styles) 1361 return self.marker_line_styles[self.imarker_line] 1362 #end def next_marker_line 1363 1364 def reset(self): 1365 self.iline = -1 1366 self.imarker = -1 1367 self.imarker_line = -1 1368 #end def reset 1369 1370 def reset_line(self): 1371 self.iline = -1 1372 #end def reset_line 1373 1374 def reset_marker(self): 1375 self.imarker = -1 1376 #end def reset_marker 1377 1378 def reset_marker_line(self): 1379 self.imarker_line = -1 1380 #end def reset_marker_line 1381#end class ColorWheel 1382color_wheel = ColorWheel() 1383 1384 1385 1386checkstats_settings = obj( 1387 verbose = False, 1388 ) 1389 1390def vlog(*args,**kwargs): 1391 if checkstats_settings.verbose: 1392 n = kwargs.get('n',0) 1393 if n==0: 1394 log(*args) 1395 else: 1396 log(*args,indent=n*' ') 1397 #end if 1398 #end if 1399#end def vlog 1400 1401 1402 1403 1404 1405# standalone definition of error function from Abramowitz & Stegun 1406# credit: http://www.johndcook.com/blog/2009/01/19/stand-alone-error-function-erf/ 1407# consider also: https://math.stackexchange.com/questions/42920/efficient-and-accurate-approximation-of-error-function 1408def erf(x): 1409 # constants 1410 a1 = 0.254829592 1411 a2 = -0.284496736 1412 a3 = 1.421413741 1413 a4 = -1.453152027 1414 a5 = 1.061405429 1415 p = 0.3275911 1416 1417 # Save the sign of x 1418 sign = 1 1419 if x < 0: 1420 sign = -1 1421 x = abs(x) 1422 1423 # A & S 7.1.26 1424 t = 1.0/(1.0 + p*x) 1425 y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*math.exp(-x*x) 1426 1427 return sign*y 1428#end def erf 1429 1430 1431# standalone inverse error function 1432# credit: https://stackoverflow.com/questions/42381244/pure-python-inverse-error-function 1433import math 1434def polevl(x, coefs, N): 1435 ans = 0 1436 power = len(coefs) - 1 1437 for coef in coefs: 1438 ans += coef * x**power 1439 power -= 1 1440 return ans 1441#end def polevl 1442 1443def p1evl(x, coefs, N): 1444 return polevl(x, [1] + coefs, N) 1445#end def p1evl 1446 1447def erfinv(z): 1448 if z < -1 or z > 1: 1449 raise ValueError("'z' must be between -1 and 1 inclusive") 1450 1451 if z == 0: 1452 return 0 1453 if z == 1: 1454 return float('inf') 1455 if z == -1: 1456 return -float('inf') 1457 1458 # From scipy special/cephes/ndrti.c 1459 def ndtri(y): 1460 # approximation for 0 <= abs(z - 0.5) <= 3/8 1461 P0 = [ 1462 -5.99633501014107895267E1, 1463 9.80010754185999661536E1, 1464 -5.66762857469070293439E1, 1465 1.39312609387279679503E1, 1466 -1.23916583867381258016E0, 1467 ] 1468 1469 Q0 = [ 1470 1.95448858338141759834E0, 1471 4.67627912898881538453E0, 1472 8.63602421390890590575E1, 1473 -2.25462687854119370527E2, 1474 2.00260212380060660359E2, 1475 -8.20372256168333339912E1, 1476 1.59056225126211695515E1, 1477 -1.18331621121330003142E0, 1478 ] 1479 1480 # Approximation for interval z = sqrt(-2 log y ) between 2 and 8 1481 # i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14. 1482 P1 = [ 1483 4.05544892305962419923E0, 1484 3.15251094599893866154E1, 1485 5.71628192246421288162E1, 1486 4.40805073893200834700E1, 1487 1.46849561928858024014E1, 1488 2.18663306850790267539E0, 1489 -1.40256079171354495875E-1, 1490 -3.50424626827848203418E-2, 1491 -8.57456785154685413611E-4, 1492 ] 1493 1494 Q1 = [ 1495 1.57799883256466749731E1, 1496 4.53907635128879210584E1, 1497 4.13172038254672030440E1, 1498 1.50425385692907503408E1, 1499 2.50464946208309415979E0, 1500 -1.42182922854787788574E-1, 1501 -3.80806407691578277194E-2, 1502 -9.33259480895457427372E-4, 1503 ] 1504 1505 # Approximation for interval z = sqrt(-2 log y ) between 8 and 64 1506 # i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890. 1507 P2 = [ 1508 3.23774891776946035970E0, 1509 6.91522889068984211695E0, 1510 3.93881025292474443415E0, 1511 1.33303460815807542389E0, 1512 2.01485389549179081538E-1, 1513 1.23716634817820021358E-2, 1514 3.01581553508235416007E-4, 1515 2.65806974686737550832E-6, 1516 6.23974539184983293730E-9, 1517 ] 1518 1519 Q2 = [ 1520 6.02427039364742014255E0, 1521 3.67983563856160859403E0, 1522 1.37702099489081330271E0, 1523 2.16236993594496635890E-1, 1524 1.34204006088543189037E-2, 1525 3.28014464682127739104E-4, 1526 2.89247864745380683936E-6, 1527 6.79019408009981274425E-9, 1528 ] 1529 1530 s2pi = 2.50662827463100050242 1531 code = 1 1532 1533 if y > (1.0 - 0.13533528323661269189): # 0.135... = exp(-2) 1534 y = 1.0 - y 1535 code = 0 1536 1537 if y > 0.13533528323661269189: 1538 y = y - 0.5 1539 y2 = y * y 1540 x = y + y * (y2 * polevl(y2, P0, 4) / p1evl(y2, Q0, 8)) 1541 x = x * s2pi 1542 return x 1543 1544 x = math.sqrt(-2.0 * math.log(y)) 1545 x0 = x - math.log(x) / x 1546 1547 z = 1.0 / x 1548 if x < 8.0: # y > exp(-32) = 1.2664165549e-14 1549 x1 = z * polevl(z, P1, 8) / p1evl(z, Q1, 8) 1550 else: 1551 x1 = z * polevl(z, P2, 8) / p1evl(z, Q2, 8) 1552 1553 x = x0 - x1 1554 if code != 0: 1555 x = -x 1556 1557 return x 1558 1559 result = ndtri((z + 1) / 2.0) / math.sqrt(2) 1560 1561 return result 1562#end def inv_erf 1563 1564 1565 1566 1567 1568# Returns failure error code to OS. 1569# Explicitly prints 'fail' after an optional message. 1570def exit_fail(msg=None): 1571 if msg!=None: 1572 print(msg) 1573 #end if 1574 print('Test status: fail') 1575 exit(1) 1576#end def exit_fail 1577 1578 1579# Returns success error code to OS. 1580# Explicitly prints 'pass' after an optional message. 1581def exit_pass(msg=None): 1582 if msg!=None: 1583 print(msg) 1584 #end if 1585 print('Test status: pass') 1586 exit(0) 1587#end def exit_pass 1588 1589 1590# Calculates the mean, variance, errorbar, and autocorrelation time 1591# for a N-d array of statistical data values. 1592# If 'exclude' is provided, the first 'exclude' values will be 1593# excluded from the analysis. 1594def simstats(x,dim=None,exclude=None): 1595 if exclude!=None: 1596 x = x[exclude:] 1597 #end if 1598 shape = x.shape 1599 ndim = len(shape) 1600 if dim==None: 1601 dim=ndim-1 1602 #end if 1603 permute = dim!=ndim-1 1604 reshape = ndim>2 1605 nblocks = shape[dim] 1606 if permute: 1607 r = list(range(ndim)) 1608 r.pop(dim) 1609 r.append(dim) 1610 permutation = tuple(r) 1611 r = list(range(ndim)) 1612 r.pop(ndim-1) 1613 r.insert(dim,ndim-1) 1614 invperm = tuple(r) 1615 x=x.transpose(permutation) 1616 shape = tuple(array(shape)[array(permutation)]) 1617 dim = ndim-1 1618 #end if 1619 if reshape: 1620 nvars = prod(shape[0:dim]) 1621 x=x.reshape(nvars,nblocks) 1622 rdim=dim 1623 dim=1 1624 else: 1625 nvars = shape[0] 1626 #end if 1627 1628 mean = x.mean(dim) 1629 var = x.var(dim) 1630 1631 N=nblocks 1632 1633 if ndim==1: 1634 i=0 1635 tempC=0.5 1636 kappa=0.0 1637 mtmp=mean 1638 if abs(var)<1e-15: 1639 kappa = 1.0 1640 else: 1641 ovar=1.0/var 1642 while (tempC>0 and i<(N-1)): 1643 kappa=kappa+2.0*tempC 1644 i=i+1 1645 #tempC=corr(i,x,mean,var) 1646 tempC = ovar/(N-i)*sum((x[0:N-i]-mtmp)*(x[i:N]-mtmp)) 1647 #end while 1648 if kappa == 0.0: 1649 kappa = 1.0 1650 #end if 1651 #end if 1652 Neff=(N+0.0)/(kappa+0.0) 1653 if (Neff == 0.0): 1654 Neff = 1.0 1655 #end if 1656 error=sqrt(var/Neff) 1657 else: 1658 error = zeros(mean.shape) 1659 kappa = zeros(mean.shape) 1660 for v in range(nvars): 1661 i=0 1662 tempC=0.5 1663 kap=0.0 1664 vtmp = var[v] 1665 mtmp = mean[v] 1666 if abs(vtmp)<1e-15: 1667 kap = 1.0 1668 else: 1669 ovar = 1.0/vtmp 1670 while (tempC>0 and i<(N-1)): 1671 i += 1 1672 kap += 2.0*tempC 1673 tempC = ovar/(N-i)*sum((x[v,0:N-i]-mtmp)*(x[v,i:N]-mtmp)) 1674 #end while 1675 if kap == 0.0: 1676 kap = 1.0 1677 #end if 1678 #end if 1679 Neff=(N+0.0)/(kap+0.0) 1680 if (Neff == 0.0): 1681 Neff = 1.0 1682 #end if 1683 kappa[v]=kap 1684 error[v]=sqrt(vtmp/Neff) 1685 #end for 1686 #end if 1687 1688 if reshape: 1689 x = x.reshape(shape) 1690 mean = mean.reshape(shape[0:rdim]) 1691 var = var.reshape(shape[0:rdim]) 1692 error = error.reshape(shape[0:rdim]) 1693 kappa = kappa.reshape(shape[0:rdim]) 1694 #end if 1695 if permute: 1696 x=x.transpose(invperm) 1697 #end if 1698 1699 return (mean,var,error,kappa) 1700#end def simstats 1701 1702 1703 1704def load_scalar_file(options,selector): 1705 output_files = options.output_files 1706 if selector=='auto': 1707 if 'dmc' in output_files: 1708 selector = 'dmc' 1709 elif 'scalar' in output_files: 1710 selector = 'scalar' 1711 else: 1712 exit_fail('could not load scalar file, no files present') 1713 #end if 1714 elif selector not in ('scalar','dmc'): 1715 exit_fail('could not load scalar file, invalid selector\ninvalid selector: {0}\nvalid options: scalar, dmc'.format(selector)) 1716 #end if 1717 if selector not in output_files: 1718 exit_fail('could not load scalar file, file is not present\nfile type requested: {0}'.format(selector)) 1719 #end if 1720 filepath = os.path.join(options.path,output_files[selector]) 1721 lt = loadtxt(filepath) 1722 if len(lt.shape)==1: 1723 lt.shape = (1,len(lt)) 1724 #end if 1725 data = lt[:,1:].transpose() 1726 fobj = open(filepath,'r') 1727 variables = fobj.readline().split()[2:] 1728 fobj.close() 1729 scalars = obj( 1730 file_type = selector, 1731 data = obj(), 1732 ) 1733 for i,var in enumerate(variables): 1734 scalars.data[var] = data[i,:] 1735 #end for 1736 return scalars 1737#end def load_scalar_file 1738 1739 1740 1741# Reads command line options. 1742def read_command_line(): 1743 try: 1744 1745 parser = OptionParser( 1746 usage='usage: %prog [options]', 1747 add_help_option=False, 1748 version='%prog 0.1' 1749 ) 1750 1751 parser.add_option('-h','--help',dest='help', 1752 action='store_true',default=False, 1753 help='Print help information and exit (default=%default).' 1754 ) 1755 parser.add_option('-p','--prefix',dest='prefix', 1756 default='qmc', 1757 help='Prefix for output files (default=%default). Can be a path including the file prefix.' 1758 ) 1759 parser.add_option('-s','--series',dest='series', 1760 default='0', 1761 help='Output series to analyze (default=%default).' 1762 ) 1763 parser.add_option('-e','--equilibration',dest='equilibration', 1764 default='0', 1765 help='Equilibration length in blocks (default=%default).' 1766 ) 1767 parser.add_option('-n','--nsigma',dest='nsigma', 1768 default='3', 1769 help='Sigma requirement for pass/fail (default=%default).' 1770 ) 1771 parser.add_option('-q','--quantity',dest='quantity', 1772 default='none', 1773 help = 'Quantity to check (required). If a non-default name for the quantity is used, pass in the quantity and name as a pair.' 1774 ) 1775 parser.add_option('-c','--npartial_sums',dest='npartial_sums', 1776 default='none', 1777 help = 'Partial sum count for the reference data (required)' 1778 ) 1779 parser.add_option('-r','--ref','--reference',dest='reference_file', 1780 default='none', 1781 help = 'Path to reference file containing full and partial sum reference information. The test fails if any full or partial sum exceeds nsigma deviation from the reference values. For cases like the density or spin density, the -f option should additionally be used (see below). For the energy density, a block by block check against relevant summed energy terms in scalar.dat or dmc.dat files is additionally made.' 1782 ) 1783 parser.add_option('-f','--fixed','--fixed_sum',dest='fixed_sum', 1784 action='store_true',default=False, 1785 help = 'Full sum of data takes on a fixed, non-stochastic value. In this case, when checking against reference data, check that each block satisfies the fixed sum condition. This is appropriate, e.g. for the electron density where the full sum of each block must equal the number of electrons. Typically the appropriate value is inferred automatically and applied by default (in other cases default=%default).' 1786 ) 1787 parser.add_option('-m','--make_ref','--make_reference',dest='make_reference', 1788 default='none', 1789 help='Used during test construction phase. Pass an integer list via -m corresponding to the number of partial sums to perform on the reference stat data followed by a series of MC step factors. The number of partial means must divide evenly into the number of stat field values for the quantity in question. The step factors relate the length of the test run (shorter) to the reference run (longer): #MC_test*factor=#MC_reference. Files containing the reference data will be produced, one for each step factor. For the partial sums, the reference sigma is increased so that the test fails with the expected probability specified by the inputted nsigma.' 1790 ) 1791 parser.add_option('-t','--plot_trace',dest='plot_trace', 1792 action='store_true',default=False, 1793 help='Plot traces of full and partial sums (default=%default).' 1794 ) 1795 parser.add_option('-v','--verbose', 1796 action='store_true',default=False, 1797 help='Print detailed information (default=%default).' 1798 ) 1799 1800 allowed_quantities = [ 1801 'density', 1802 'spindensity', 1803 'energydensity', 1804 '1rdm', 1805 '1redm', 1806 'momentum', 1807 ] 1808 1809 opt,files_in = parser.parse_args() 1810 options = obj() 1811 options.transfer_from(opt.__dict__) 1812 1813 if options.help: 1814 print('\n'+parser.format_help().strip()) 1815 print('\n\nExample usage:') 1816 print('\n Making reference data to create a test:') 1817 print(" check_stats.py -p qmc -s 0 -q spindensity -e 10 -c 8 -v -m '0 10 100'") 1818 print('\n Using reference data to perform a test:') 1819 print(' check_stats.py -p qmc -s 0 -q spindensity -e 10 -c 8 -n 3 -r qmc.s000.stat_ref_spindensity_10.dat') 1820 print() 1821 exit() 1822 #end if 1823 1824 if len(files_in)>0: 1825 exit_fail('check_stats does not accept file as input, only command line arguments\nfiles provided: {0}'.format(files_in)) 1826 #end if 1827 1828 checkstats_settings.verbose = options.verbose 1829 1830 vlog('\nreading command line inputs') 1831 1832 options.series = int(options.series) 1833 options.equilibration = int(options.equilibration) 1834 options.nsigma = float(options.nsigma) 1835 options.path,options.prefix = os.path.split(options.prefix) 1836 1837 if options.plot_trace and not can_plot: 1838 vlog('trace plots requested, but plotting libraries are not available\ndisabling plots',n=1) 1839 options.plot_trace = False 1840 #end if 1841 1842 if options.path=='': 1843 options.path = './' 1844 #end if 1845 1846 options.qlabel = None 1847 if ' ' in options.quantity or ',' in options.quantity: 1848 qlist = options.quantity.strip('"').strip("'").replace(',',' ').split() 1849 if len(qlist)!=2: 1850 exit_fail('quantity can accept only one or two values\nyou provided {0}: {1}'.format(len(qlist),qlist)) 1851 #end if 1852 options.quantity,options.qlabel = qlist 1853 #end if 1854 if options.qlabel is None: 1855 default_label = obj({ 1856 'density' : 'Density' , 1857 'spindensity' : 'SpinDensity' , 1858 'energydensity' : 'EnergyDensity' , 1859 '1rdm' : 'DensityMatrices', 1860 '1redm' : 'DensityMatrices', 1861 'momentum' : 'nofk' , 1862 }) 1863 options.qlabel = default_label[options.quantity] 1864 #end if 1865 if options.quantity=='none': 1866 exit_fail('must provide quantity') 1867 elif options.quantity not in allowed_quantities: 1868 exit_fail('unrecognized quantity provided\nallowed quantities: {0}\nquantity provided: {1}'.format(allowed_quantities,options.quantity)) 1869 #end if 1870 1871 if options.npartial_sums=='none': 1872 exit_fail('-c option is required') 1873 #end if 1874 options.npartial_sums = int(options.npartial_sums) 1875 1876 if options.reference_file!='none': 1877 if not os.path.exists(options.reference_file): 1878 exit_fail('reference file does not exist\nreference file provided: {0}'.format(options.reference_file)) 1879 #end if 1880 options.make_reference = False 1881 elif options.make_reference!='none': 1882 try: 1883 mr = array(options.make_reference.split(),dtype=int) 1884 except: 1885 exit_fail('make_reference must be a list of integers\nyou provided: {0}'.format(options.make_reference)) 1886 #end try 1887 if len(mr)<1: 1888 exit_fail('make_reference must contain at least one MC length factor') 1889 #end if 1890 options.mc_factors = mr 1891 options.make_reference = True 1892 else: 1893 exit_fail('must provide either reference_file or make_reference') 1894 #end if 1895 1896 fixed_sum_quants = set(['density','spindensity','energydensity']) 1897 if options.quantity in fixed_sum_quants: 1898 options.fixed_sum = True 1899 #end if 1900 1901 vlog('inputted options:\n'+str(options),n=1) 1902 1903 except Exception as e: 1904 exit_fail('error during command line read:\n'+str(e)) 1905 #end try 1906 1907 return options 1908#end def read_command_line 1909 1910 1911 1912 1913def process_stat_file(options): 1914 vlog('processing stat.h5 file') 1915 1916 values = obj() 1917 1918 try: 1919 # find all output files matching prefix 1920 vlog('searching for qmcpack output files',n=1) 1921 vlog('search path:\n '+options.path,n=2) 1922 prefix = options.prefix+'.s'+str(options.series).zfill(3) 1923 files = os.listdir(options.path) 1924 output_files = obj() 1925 for file in files: 1926 if file.startswith(prefix): 1927 if file.endswith('.stat.h5'): 1928 output_files.stat = file 1929 elif file.endswith('.scalar.dat'): 1930 output_files.scalar = file 1931 elif file.endswith('.dmc.dat'): 1932 output_files.dmc = file 1933 #end if 1934 #end if 1935 #end for 1936 options.output_files = output_files 1937 vlog('files found:\n'+str(output_files).rstrip(),n=2) 1938 if 'stat' not in output_files: 1939 exit_fail('stat.h5 file matching prefix {0} was not found\nsearch path: {1}'.format(prefix,options.path)) 1940 #end if 1941 1942 # read data from the stat file 1943 vlog('opening stat.h5 file',n=1) 1944 stat = read_hdf(os.path.join(options.path,output_files.stat),view=True) 1945 vlog('file contents:\n'+repr(stat).rstrip(),n=2) 1946 vlog('extracting {0} data'.format(options.quantity),n=1) 1947 vlog('searching for {0} with label {1}'.format(options.quantity,options.qlabel),n=2) 1948 if options.qlabel in stat: 1949 qstat = stat[options.qlabel] 1950 vlog('{0} data contents:\n{1}'.format(options.quantity,repr(qstat).rstrip()),n=2) 1951 else: 1952 exit_fail('could not find {0} data with label {1}'.format(options.quantity,options.qlabel)) 1953 #end if 1954 quantity_paths = obj({ 1955 'density' : obj(tot='value'), 1956 'spindensity' : obj(u='u/value', 1957 d='d/value'), 1958 '1rdm' : obj(u='number_matrix/u/value', 1959 d='number_matrix/d/value'), 1960 '1redm' : obj(u='energy_matrix/u/value', 1961 d='energy_matrix/d/value'), 1962 'energydensity' : obj(W=('spacegrid1/value',0,3), 1963 T=('spacegrid1/value',1,3), 1964 V=('spacegrid1/value',2,3)), 1965 'momentum' : obj(tot='value'), 1966 }) 1967 qpaths = quantity_paths[options.quantity] 1968 vlog('search paths:\n{0}'.format(str(qpaths).rstrip()),n=2) 1969 qdata = obj() 1970 dfull = None 1971 for dname,dpath in qpaths.items(): 1972 packed = isinstance(dpath,tuple) 1973 if packed: 1974 dpath,dindex,dcount = dpath 1975 #end if 1976 if not qstat.path_exists(dpath): 1977 exit_fail('{0} data not found in file {1}\npath searched: {2}'.format(options.quantity,output_files.stat,dpath)) 1978 #end if 1979 if not packed: 1980 d = array(qstat.get_path(dpath),dtype=float) 1981 else: 1982 if dfull is None: 1983 dfull = array(qstat.get_path(dpath),dtype=float) 1984 dfull.shape = dfull.shape[0],dfull.shape[1]//dcount,dcount 1985 #end if 1986 d = dfull[:,:,dindex] 1987 d.shape = dfull.shape[0],dfull.shape[1] 1988 #end if 1989 qdata[dname] = d 1990 vlog('{0} data found with shape {1}'.format(dname,d.shape),n=2) 1991 if len(d.shape)>2: 1992 d.shape = d.shape[0],d.size//d.shape[0] 1993 vlog('reshaped {0} data to {1}'.format(dname,d.shape),n=2) 1994 #end if 1995 options.nblocks = d.shape[0] 1996 #end for 1997 1998 # process the data, taking full and partial sums 1999 vlog('processing {0} data'.format(options.quantity),n=1) 2000 for dname,d in qdata.items(): 2001 vlog('processing {0} data'.format(dname),n=2) 2002 if d.shape[1]%options.npartial_sums!=0: 2003 exit_fail('cannot make partial sums\nnumber of requested partial sums does not divide evenly into the number of values available\nrequested partial sums: {0}\nnumber of values present: {1}\nnvalue/npartial_sums: {2}'.format(options.npartial_sums,d.shape[1],float(d.shape[1])/options.npartial_sums)) 2004 #end if 2005 data = obj() 2006 data.full_sum = d.sum(1) 2007 vlog('full sum data shape: {0}'.format(data.full_sum.shape),n=3) 2008 data.partial_sums = zeros((d.shape[0],options.npartial_sums)) 2009 psize = d.shape[1]//options.npartial_sums 2010 for p in range(options.npartial_sums): 2011 data.partial_sums[:,p] = d[:,p*psize:(p+1)*psize].sum(1) 2012 #end for 2013 vlog('partial sum data shape: {0}'.format(data.partial_sums.shape),n=3) 2014 fmean,var,ferror,kappa = simstats(data.full_sum,exclude=options.equilibration) 2015 vlog('full sum mean : {0}'.format(fmean),n=3) 2016 vlog('full sum error: {0}'.format(ferror),n=3) 2017 pmean,var,perror,kappa = simstats(data.partial_sums,dim=0,exclude=options.equilibration) 2018 vlog('partial sum mean : {0}'.format(pmean),n=3) 2019 vlog('partial sum error: {0}'.format(perror),n=3) 2020 values[dname] = obj( 2021 full_mean = fmean, 2022 full_error = ferror, 2023 partial_mean = pmean, 2024 partial_error = perror, 2025 data = data, 2026 ) 2027 #end for 2028 2029 # check that all values have been processed 2030 missing = set(qpaths.keys())-set(values.keys()) 2031 if len(missing)>0: 2032 exit_fail('some values not processed\nvalues missing: {0}'.format(sorted(missing))) 2033 #end if 2034 2035 # plot quantity traces, if requested 2036 if options.plot_trace: 2037 vlog('creating trace plots of full and partial sums',n=1) 2038 for dname,dvalues in values.items(): 2039 label = options.quantity 2040 if len(values)>1: 2041 label+=' '+dname 2042 #end if 2043 data = dvalues.data 2044 figure() 2045 plot(data.full_sum) 2046 title('Trace of {0} full sum'.format(label)) 2047 xlabel('Block index') 2048 figure() 2049 plot(data.partial_sums) 2050 title('Trace of {0} partial sums'.format(label)) 2051 xlabel('Block index') 2052 #end for 2053 show() 2054 #end if 2055 except Exception as e: 2056 exit_fail('error during stat file processing:\n'+str(e)) 2057 #end try 2058 2059 return values 2060#end def process_stat_file 2061 2062 2063 2064def make_reference_files(options,values): 2065 vlog('\nmaking reference files') 2066 2067 # create a reference file for each Monte Carlo sample factor 2068 for mcfac in options.mc_factors: 2069 errfac = sqrt(1.0+mcfac) 2070 filename = '{0}.s{1}.stat_ref_{2}_{3}.dat'.format(options.prefix,str(options.series).zfill(3),options.quantity,mcfac) 2071 filepath = os.path.join(options.path,filename) 2072 vlog('writing reference file for {0}x shorter test runs'.format(mcfac),n=1) 2073 vlog('reference file location: '+filepath,n=2) 2074 f = open(filepath,'w') 2075 # write descriptive header line 2076 line = '# ' 2077 for dname in sorted(values.keys()): 2078 line += ' {0:<16} {1:<16}'.format(dname,dname+'_err') 2079 #end for 2080 f.write(line+'\n') 2081 # write means and errors of full sum 2082 line = '' 2083 for dname in sorted(values.keys()): 2084 dvalues = values[dname] 2085 fmean = dvalues.full_mean 2086 ferror = dvalues.full_error 2087 line += ' {0: 16.8e} {1: 16.8e}'.format(fmean,errfac*ferror) 2088 #end for 2089 f.write(line+'\n') 2090 # write means and errors of partial sums 2091 for p in range(options.npartial_sums): 2092 line = '' 2093 for dname in sorted(values.keys()): 2094 dvalues = values[dname] 2095 pmean = dvalues.partial_mean 2096 perror = dvalues.partial_error 2097 line += ' {0: 16.8e} {1: 16.8e}'.format(pmean[p],errfac*perror[p]) 2098 #end for 2099 f.write(line+'\n') 2100 #end for 2101 f.close() 2102 #end for 2103 2104 # create a trace file containing full and partial sum data per block 2105 filename = '{0}.s{1}.stat_trace_{2}.dat'.format(options.prefix,str(options.series).zfill(3),options.quantity) 2106 filepath = os.path.join(options.path,filename) 2107 vlog('writing trace file containing full and partial sums per block',n=1) 2108 vlog('trace file location: '+filepath,n=2) 2109 f = open(filepath,'w') 2110 # write descriptive header line 2111 line = '# index ' 2112 for dname in sorted(values.keys()): 2113 line += ' {0:<16}'.format(dname+'_full') 2114 for p in range(options.npartial_sums): 2115 line += ' {0:<16}'.format(dname+'_partial_'+str(p)) 2116 #end for 2117 #end for 2118 f.write(line+'\n') 2119 # write full and partial sum data per block 2120 for b in range(options.nblocks): 2121 line = ' {0:>6}'.format(b) 2122 for dname in sorted(values.keys()): 2123 dvalues = values[dname].data 2124 fsum = dvalues.full_sum 2125 psums = dvalues.partial_sums[b] 2126 line += ' {0: 16.8e}'.format(fsum[b]) 2127 for psum in psums: 2128 line += ' {0: 16.8e}'.format(psum) 2129 #end for 2130 #end for 2131 f.write(line+'\n') 2132 #end for 2133 f.close() 2134 vlog('\n') 2135#end def make_reference_files 2136 2137 2138 2139# Checks computed values from stat.h5 files against specified reference values. 2140passfail = {True:'pass',False:'fail'} 2141def check_values(options,values): 2142 vlog('\nchecking against reference values') 2143 2144 success = True 2145 msg = '' 2146 2147 try: 2148 msg += '\nTests for series {0} quantity "{1}"\n'.format(options.series,options.quantity) 2149 2150 # find nsigma for each partial sum 2151 # overall probability of partial sum failure is according to original nsigma 2152 vlog('adjusting nsigma to account for partial sum count',n=1) 2153 x = longdouble(options.nsigma/sqrt(2.)) 2154 N = options.npartial_sums 2155 nsigma_partial = sqrt(2.)*erfinv(erf(x)**(1./N)) 2156 vlog('overall full/partial test nsigma: {0}'.format(options.nsigma),n=2) 2157 vlog('adjusted per partial sum nsigma : {0}'.format(nsigma_partial),n=2) 2158 2159 # read in the reference file 2160 vlog('reading reference file',n=1) 2161 vlog('reference file location: {0}'.format(options.reference_file),n=2) 2162 f = open(options.reference_file,'r') 2163 dnames = f.readline().split()[1::2] 2164 vlog('sub-quantities found: {0}'.format(dnames),n=2) 2165 if set(dnames)!=set(values.keys()): 2166 missing = set(values.keys())-set(dnames) 2167 extra = set(dnames)-set(values.keys()) 2168 if missing>0: 2169 exit_fail('some sub-quantities are missing\npresent in test files: {0}\npresent in reference files: {1}\nmissing: {2}'.format(sorted(values.keys()),sorted(dnames),sorted(missing))) 2170 elif extra>0: 2171 exit_fail('some sub-quantities are extra\npresent in test files: {0}\npresent in reference files: {1}\nextra: {2}'.format(sorted(values.keys()),sorted(dnames),sorted(extra))) 2172 else: 2173 exit_fail('developer error, this point should be impossible to reach') 2174 #end if 2175 #end if 2176 ref = array(f.read().split(),dtype=float) 2177 ref.shape = len(ref)//(2*len(dnames)),2*len(dnames) 2178 full = ref[0,:].ravel() 2179 partial = ref[1:,:].T 2180 if len(ref)-1!=options.npartial_sums: 2181 exit_fail('test and reference partial sum counts do not match\ntest partial sum count: {0}\nreference partial sum count: {1}'.format(options.npartial_sums,len(ref)-1)) 2182 #end if 2183 vlog('partial sum count found: {0}'.format(len(ref)-1),n=2) 2184 ref_values = obj() 2185 for dname in dnames: 2186 ref_values[dname] = obj() 2187 #end for 2188 n=0 2189 for dname in dnames: 2190 ref_values[dname].set( 2191 full_mean = full[n], 2192 full_error = full[n+1], 2193 partial_mean = partial[n,:].ravel(), 2194 partial_error = partial[n+1,:].ravel(), 2195 ) 2196 n+=2 2197 #end for 2198 npartial = len(ref)-1 2199 f.close() 2200 dnames = sorted(dnames) 2201 vlog('reference file read successfully',n=2) 2202 2203 # for cases with fixed full sum, check the per block sum 2204 if options.fixed_sum: 2205 vlog('checking per block fixed sums',n=1) 2206 msg+='\n Fixed sum per block tests:\n' 2207 dnames_fixed = dnames 2208 if options.quantity=='energydensity': 2209 dnames_fixed = ['W'] 2210 #end if 2211 fixed_sum_success = True 2212 ftol = 1e-8 2213 eq = options.equilibration 2214 for dname in dnames_fixed: 2215 ref_vals = ref_values[dname] 2216 ref_mean = ref_vals.full_mean 2217 ref_error = ref_vals.full_error 2218 if abs(ref_error/ref_mean)>ftol: 2219 exit_fail('reference fixed sum is not fixed as asserted\ncannot check per block fixed sums\nplease check reference data') 2220 #end if 2221 test_vals = values[dname].data.full_sum 2222 for i,v in enumerate(test_vals[eq:]): 2223 if abs((v-ref_mean)/ref_mean)>ftol: 2224 fixed_sum_success = False 2225 msg += ' {0} {1} {2}!={3}\n'.format(dname,i,v,ref_mean) 2226 #end if 2227 #end for 2228 #end for 2229 if fixed_sum_success: 2230 fmsg = 'all per block sums match the reference' 2231 else: 2232 fmsg = 'some per block sums do not match the reference' 2233 #end if 2234 vlog(fmsg,n=2) 2235 msg += ' '+fmsg+'\n' 2236 msg += ' status of this test: {0}\n'.format(passfail[fixed_sum_success]) 2237 success &= fixed_sum_success 2238 #end if 2239 2240 # for the energy density, check per block sums against the scalar file 2241 if options.quantity=='energydensity': 2242 vlog('checking energy density terms per block',n=1) 2243 msg+='\n Energy density sums vs. scalar file per block tests:\n' 2244 scalars = load_scalar_file(options,'auto') 2245 ed_success = True 2246 ftol = 1e-8 2247 ed_values = obj( 2248 T = values.T.data.full_sum, 2249 V = values.V.data.full_sum, 2250 ) 2251 ed_values.E = ed_values.T + ed_values.V 2252 if scalars.file_type=='scalar': 2253 comparisons = obj( 2254 E='LocalEnergy', 2255 T='Kinetic', 2256 V='LocalPotential', 2257 ) 2258 elif scalars.file_type=='dmc': 2259 comparisons = obj(E='LocalEnergy') 2260 else: 2261 exit_fail('unrecognized scalar file type: {0}'.format(scalars.file_type)) 2262 #end if 2263 for k in sorted(comparisons.keys()): 2264 ed_vals = ed_values[k] 2265 sc_vals = scalars.data[comparisons[k]] 2266 if scalars.file_type=='dmc': 2267 if len(sc_vals)%len(ed_vals)==0 and len(sc_vals)>len(ed_vals): 2268 steps = len(sc_vals)//len(ed_vals) 2269 sc_vals.shape = len(sc_vals)//steps,steps 2270 sc_vals = sc_vals.mean(1) 2271 #end if 2272 #end if 2273 if len(ed_vals)!=len(sc_vals): 2274 exit_fail('energy density per block test cannot be completed\nnumber of energy density and scalar blocks do not match\nenergy density blocks: {0}\nscalar file blocks: {1}'.format(len(ed_vals),len(sc_vals))) 2275 #end if 2276 for i,(edv,scv) in enumerate(zip(ed_vals,sc_vals)): 2277 if abs((edv-scv)/scv)>ftol: 2278 ed_success = False 2279 msg += ' {0} {1} {2}!={3}\n'.format(k,i,edv,scv) 2280 #end if 2281 #end for 2282 #end for 2283 if ed_success: 2284 fmsg = 'all per block sums match the scalar file' 2285 else: 2286 fmsg = 'some per block sums do not match the scalar file' 2287 #end if 2288 vlog(fmsg,n=2) 2289 msg += ' '+fmsg+'\n' 2290 msg += ' status of this test: {0}\n'.format(passfail[ed_success]) 2291 success &= ed_success 2292 #end if 2293 2294 2295 # function used immediately below to test a mean value vs reference 2296 def check_mean(label,mean_comp,error_comp,mean_ref,error_ref,nsigma): 2297 msg='\n Testing quantity: {0}\n'.format(label) 2298 2299 # ensure error_ref is large enough for non-statistical quantities 2300 ctol = 1e-12 2301 if abs(error_ref/mean_ref)<ctol: 2302 error_ref = ctol*mean_ref 2303 #end if 2304 2305 quant_success = abs(mean_comp-mean_ref) <= nsigma*error_ref 2306 2307 delta = mean_comp-mean_ref 2308 delta_err = sqrt(error_comp**2+error_ref**2) 2309 2310 msg+=' reference mean value : {0: 12.8f}\n'.format(mean_ref) 2311 msg+=' reference error bar : {0: 12.8f}\n'.format(error_ref) 2312 msg+=' computed mean value : {0: 12.8f}\n'.format(mean_comp) 2313 msg+=' computed error bar : {0: 12.8f}\n'.format(error_comp) 2314 msg+=' pass tolerance : {0: 12.8f} ({1: 12.8f} sigma)\n'.format(nsigma*error_ref,nsigma) 2315 if error_ref > 0.0: 2316 msg+=' deviation from reference : {0: 12.8f} ({1: 12.8f} sigma)\n'.format(delta,delta/error_ref) 2317 #end if 2318 msg+=' error bar of deviation : {0: 12.8f}\n'.format(delta_err) 2319 if error_ref > 0.0: 2320 msg+=' significance probability : {0: 12.8f} (gaussian statistics)\n'.format(erf(abs(delta/error_ref)/math.sqrt(2.0))) 2321 #end if 2322 msg+=' status of this test : {0}\n'.format(passfail[quant_success]) 2323 2324 return quant_success,msg 2325 #end def check_mean 2326 2327 2328 # check full and partial sums vs the reference 2329 vlog('checking full and partial sums',n=1) 2330 for dname in dnames: 2331 vals = values[dname] 2332 ref_vals = ref_values[dname] 2333 # check full sum 2334 vlog('checking full sum mean for "{0}"'.format(dname),n=2) 2335 qsuccess,qmsg = check_mean( 2336 label = '{0} full sum'.format(dname), 2337 mean_comp = vals.full_mean, 2338 error_comp = vals.full_error, 2339 mean_ref = ref_vals.full_mean, 2340 error_ref = ref_vals.full_error, 2341 nsigma = options.nsigma, 2342 ) 2343 vlog('status for full sum: {0}'.format(passfail[qsuccess]),n=3) 2344 msg += qmsg 2345 success &= qsuccess 2346 # check partial sums 2347 vlog('checking partial sum means for "{0}"'.format(dname),n=2) 2348 for p in range(len(ref_vals.partial_mean)): 2349 qsuccess,qmsg = check_mean( 2350 label = '{0} partial sum {1}'.format(dname,p), 2351 mean_comp = vals.partial_mean[p], 2352 error_comp = vals.partial_error[p], 2353 mean_ref = ref_vals.partial_mean[p], 2354 error_ref = ref_vals.partial_error[p], 2355 nsigma = nsigma_partial, 2356 ) 2357 vlog('status for partial sum {0}: {1}'.format(p,passfail[qsuccess]),n=3) 2358 msg += qmsg 2359 success &= qsuccess 2360 #end for 2361 #end for 2362 except Exception as e: 2363 exit_fail('error during value check:\n'+str(e)) 2364 #end try 2365 2366 return success,msg 2367#end def check_values 2368 2369 2370 2371 2372# Main execution 2373if __name__=='__main__': 2374 # Read and interpret command line options. 2375 options = read_command_line() 2376 2377 # Compute means of desired quantities from a stat.h5 file. 2378 values = process_stat_file(options) 2379 2380 if options.make_reference: 2381 # Make reference files to create tests from. 2382 make_reference_files(options,values) 2383 exit() 2384 else: 2385 # Check computed means against reference solutions. 2386 success,msg = check_values(options,values) 2387 2388 # Pass success/failure exit codes and strings to the OS. 2389 if success: 2390 exit_pass(msg) 2391 else: 2392 exit_fail(msg) 2393 #end if 2394 #end if 2395#end if 2396