Changeset 903 for trunk/GSASIImapvars.py
 Timestamp:
 May 13, 2013 3:18:21 PM (8 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIImapvars.py
r762 r903 7 7 # $Id$ 8 8 ########### SVN repository information ################### 9 """Module that implements algebraic contraints, parameter redefinition 9 """ 10 *GSASIImapvars: Parameter constraints* 11 ====================================== 12 13 Module to implements algebraic contraints, parameter redefinition 10 14 and parameter simplification contraints. 11 15 12 16 Parameter redefinition (new vars) is done by creating one or more relationships 13 17 between a set of parameters 18 19 :: 20 14 21 Mx1 * Px + My1 * Py +... 15 22 Mx2 * Px + Mz2 * Pz + ... 23 16 24 where Pj is a parameter name and Mjk is a constant. 17 25 18 26 Constant constraint Relations can also be supplied in the form of an equation: 27 28 :: 29 19 30 nx1 * Px + ny1 * Py +... = C1 31 20 32 where Cn is a constant. These equations define an algebraic 21 33 constrant. … … 23 35 Parameters can also be "fixed" (held), which prevents them from being refined. 24 36 25 All of the above three cases is supplied asinput using routines37 All of the above three cases are input using routines 26 38 GroupConstraints and GenerateConstraints. The input consists of a list of 27 39 relationship dictionaries: 40 41 .. codeblock:: python 42 28 43 constrDict = [ 29 44 {'0:12:Scale': 2.0, '0:14:Scale': 4.0, '0:13:Scale': 3.0, '0:0:Scale': 0.5}, 30 45 {'2::C(10,6,1)': 1.0, '1::C(10,6,1)': 1.0}, 31 46 {'0::A0': 0.0}] 32 33 47 fixedList = ['5.0', None, '0'] 34 48 … … 41 55 (independent) parameter is equated to several (now dependent) parameters. In 42 56 algebraic form this is: 57 58 :: 59 43 60 P0 = M1 * P1 = M2 * P2 = ... 61 44 62 Thus parameters P0, P1 and P2,... are linearly equivalent. Routine StoreEquivalence is 45 63 used to specify these equivalences. … … 68 86 varied (appear in varyList) or none can be varied. This is checked in GenerateConstraints 69 87 (as well as for generated relationships in SetVaryFlags). 88 70 89 * If all parameters in a parameter redefinition (new var) relationship are varied, the 71 90 parameter assigned to this expression (::constr:n, see paramPrefix) newly generated 72 91 parameter is varied. Note that any generated "missing" relations are not varied. Only 73 92 the input relations are varied. 93 74 94 * If all parameters in a fixed constraint equation are varied, the generated "missing" 75 95 relations in the group are all varied. This provides the NC degrees of freedom. 76 96 77 External Routines: 78 To define a set of constrained and unconstrained relations, one 79 calls InitVars, GroupConstraints and GenerateConstraints. It 80 is best to supply all relations/equations in one call to these 81 routines so that grouping is done correctly. 82 83 To implement parameter redefinition, one calls 84 StoreEquivalence. This should be called for every set of 85 equivalence relationships. There is no harm in using 86 StoreEquivalence with the same independent variable: 97 *External Routines* 98  99 100 To define a set of constrained and unconstrained relations, one 101 defines a list of dictionary defining constraint parameters and their 102 values, a list of fixed values for each constraint and a list of 103 parameters to be varied. In addition, one uses 104 :func:`StoreEquivalence` to define parameters that are equivalent. One 105 can then use :func:`CheckConstraints` to check that the input is 106 internally consistent and finally :func:`GroupConstraints` and 107 :func:`GenerateConstraints` to generate the internally used 108 tables. Routines :func:`Map2Dict` is used to initialize the parameter 109 dictionary and :func:`Dict2Map`, :func:`Dict2Deriv`, and 110 :func:`ComputeDepESD` are used to apply constraints. Routine 111 :func:`VarRemapShow` is used to print out the constraint information, 112 as stored by :func:`GenerateConstraints`. 113 114 :func:`InitVars` 115 This is optionally used to clear out all defined previously defined constraint information 116 117 :func:`StoreEquivalence` 118 To implement parameter redefinition, one calls StoreEquivalence. This should be called for every set of 119 equivalence relationships. There is no harm in using StoreEquivalence with the same independent variable: 120 121 .. codeblock:: python 122 87 123 StoreEquivalence('x',('y',)) 88 124 StoreEquivalence('x',('z',)) 89 (though StoreEquivalence('x',('y','z')) will run more 90 efficiently) but mixing independent and dependent variables is 91 problematic. This is not allowed: 125 126 or equivalently 127 128 .. codeblock:: python 129 130 StoreEquivalence('x',('y','z')) 131 132 The latter will run more efficiently. Note that mixing independent and dependent variables is 133 problematic. This is not allowed: 134 135 .. codeblock:: python 136 92 137 StoreEquivalence('x',('y',)) 93 138 StoreEquivalence('y',('z',)) 94 Use StoreEquivalence before calling GenerateConstraints or 95 CheckConstraints 96 139 140 Use StoreEquivalence before calling GenerateConstraints or CheckConstraints 141 142 :func:`CheckConstraints` 97 143 To check that input in internally consistent, use CheckConstraints 98 144 99 To show a summary of the parameter remapping, one calls 100 VarRemapShow 101 145 :func:`Map2Dict` 102 146 To determine values for the parameters created in this module, one 103 calls Map2Dict. This will not apply contraints. 104 147 calls Map2Dict. This will not apply contraints. 148 149 :func:`Dict2Map` 105 150 To take values from the new independent parameters and constraints, 106 one calls Dict2Map. This will apply contraints. 107 151 one calls Dict2Map. This will apply contraints. 152 153 :func:`Dict2Deriv` 108 154 Use Dict2Deriv to determine derivatives on independent parameters 109 from those on dependent ones 110 155 from those on dependent ones 156 157 :func:`ComputeDepESD` 111 158 Use ComputeDepESD to compute uncertainties on dependent variables 112 159 113 Global Variables: 114 dependentParmList: contains a list by group of lists of 115 parameters used in the group. Note that parameters listed in 116 dependentParmList should not be refined as they will not affect 117 the model 118 indParmList: a list of groups of Independent parameters defined in 160 :func:`VarRemapShow` 161 To show a summary of the parameter remapping, one calls VarRemapShow 162 163 *Global Variables* 164  165 166 dependentParmList: 167 contains a list by group of lists of 168 parameters used in the group. Note that parameters listed in 169 dependentParmList should not be refined as they will not affect 170 the model 171 172 indParmList: 173 a list of groups of Independent parameters defined in 119 174 each group. This contains both parameters used in parameter 120 175 redefinitions as well as names of generated new parameters. 121 176 122 fixedVarList: a list of variables that have been 'fixed' 177 fixedVarList: 178 a list of variables that have been 'fixed' 123 179 by defining them as equal to a constant (::var: = 0). Note that 124 180 the constant value is ignored at present. These variables are 125 181 later removed from varyList which prevents them from being refined. 126 182 Unlikely to be used externally. 127 arrayList: a list by group of relationship matrices to relate 183 184 arrayList: 185 a list by group of relationship matrices to relate 128 186 parameters in dependentParmList to those in indParmList. Unlikely 129 187 to be used externally. 130 invarrayList: a list by group of relationship matrices to relate 188 189 invarrayList: 190 a list by group of relationship matrices to relate 131 191 parameters in indParmList to those in dependentParmList. Unlikely 132 192 to be used externally. 133 fixedDict: a dictionary containing the fixed values corresponding 193 194 fixedDict: 195 a dictionary containing the fixed values corresponding 134 196 to parameter equations. The dict key is an ascii string, but the 135 197 dict value is a float. Unlikely to be used externally. 198 199 *Routines* 200  201 202 Note that parameter names in GSASII are strings of form ``<ph>:<hst>:<nam>`` 203 136 204 """ 137 205 … … 167 235 def GroupConstraints(constrDict): 168 236 """divide the constraints into groups that share no parameters. 169 Input 170 constrDict: a list of dicts defining relationships/constraints 237 238 :param dict constrDict: a list of dicts defining relationships/constraints 239 240 :: 241 171 242 constrDict = [{<constr1>}, {<constr2>}, ...] 172 where {<constr1>} is {'param1': mult1, 'param2': mult2,...} 173 Returns two lists of lists: 174 a list of relationship list indices for each group pointing to 175 elements in constrDict and 176 a list containing the parameters used in each group 243 244 where {<constr1>} is {'param1': mult1, 'param2': mult2,...} 245 246 :returns: two lists of lists: 247 248 * a list of grouped contraints where each constraint grouped containts a list of indices for constraint constrDict entries 249 * a list containing lists of parameter names contained in each group 250 177 251 """ 178 252 assignedlist = [] # relationships that have been used … … 205 279 constraints and checks for inconsistencies such as conflicts in 206 280 parameter/variable definitions and or inconsistently varied parameters. 207 Input: see GenerateConstraints 208 Output: returns two strings 209 the first lists conflicts internal to the specified constraints 210 the second lists conflicts where the varyList specifies some 281 282 :param list varyList: a list of parameters names that will be varied 283 284 :param dict constrDict: a list of dicts defining relationships/constraints (as defined in :func:`GroupConstraints`) 285 286 :param list fixedList: a list of values specifying a fixed value for each dict in constrDict. Values are 287 either strings that can be converted to floats or None if the constraint defines a new parameter rather 288 than a constant. 289 290 :returns: two strings: 291 292 * the first lists conflicts internal to the specified constraints 293 * the second lists conflicts where the varyList specifies some 211 294 parameters in a constraint, but not all 295 212 296 If there are no errors, both strings will be empty 213 297 ''' … … 267 351 errmsg += str(mv) + " => " + s + '\n' 268 352 269 #print 'indepVarList',indepVarList270 #print 'depVarList',depVarList271 353 # check for errors: 272 354 inboth = set(indepVarList).intersection(set(depVarList)) … … 313 395 if var in fixVlist: 314 396 errmsg += '\nParameter '+var+" is Fixed and used in a constraint:\n\t" 315 errmsg += FormatConstraint(constrDict[rel],fixedList[rel])+"\n" 316 # if var in equivVarList: 317 # errmsg += '\nParameter '+var+" is Equivalenced and used in a constraint:\n\t" 318 # errmsg += FormatConstraint(constrDict[rel],fixedList[rel])+"\n" 397 errmsg += _FormatConstraint(constrDict[rel],fixedList[rel])+"\n" 319 398 if varied > 0 and varied != len(constrDict[rel]): 320 399 warnmsg += "\nNot all variables refined in constraint:\n\t" 321 warnmsg += FormatConstraint(constrDict[rel],fixedList[rel])400 warnmsg += _FormatConstraint(constrDict[rel],fixedList[rel]) 322 401 warnmsg += '\nNot refined: ' + notvaried + '\n' 323 if errmsg or warnmsg: return errmsg,warnmsg 402 if errmsg or warnmsg: 403 return errmsg,warnmsg 324 404 325 405 # now look for process each group and create the relations that are needed to form 326 406 # nonsingular square matrix 327 407 for group,varlist in zip(groups,parmlist): 328 if len(varlist) == 1: continue 329 try: 330 arr = MakeArray(constrDict, group, varlist) 408 if len(varlist) == 1: continue # a constraint group with a single variable can be ignored 409 if len(varlist) < len(group): # too many relationships  no can do 410 errmsg += "\nOverconstrained input. " 411 errmsg += "There are more constraints " + str(len(group)) 412 errmsg += "\n\tthan variables " + str(len(varlist)) + "\n" 413 for rel in group: 414 errmsg += _FormatConstraint(constrDict[rel],fixedList[rel]) 415 errmsg += "\n" 416 continue 417 try: 418 multarr = _FillArray(group,constrDict,varlist) 419 _RowEchelon(len(group),multarr,varlist) 331 420 except: 332 errmsg += "\nOverconstrained input. " 333 errmsg += "There are more constraints" + str(len(group)) 334 errmsg += "\n\tthan variables" + str(len(varlist)) + "\n" 421 errmsg += "\nSingular input. " 422 errmsg += "There are internal inconsistencies in these constraints\n" 335 423 for rel in group: 336 errmsg += FormatConstraint(constrDict[rel],fixedList[rel])424 errmsg += _FormatConstraint(constrDict[rel],fixedList[rel]) 337 425 errmsg += "\n" 426 continue 338 427 try: 339 constrArr = FillInMissingRelations(arr,len(group)) 340 except Exception,msg: 341 if msg.message.startswith('VectorProduct'): 342 errmsg += "\nSingular input. " 343 errmsg += "This set of constraints is not linearly independent:\n\n" 344 else: 345 errmsg += "\nInconsistent input. " 346 errmsg += "Cound not generate a set of nonsingular constraints" 347 errmsg += "\n\tfor this set of input:\n\n" 428 multarr = _FillArray(group,constrDict,varlist,FillDiagonals=True) 429 GramSchmidtOrtho(multarr,len(group)) 430 except: 431 errmsg += "\nUnexpected singularity with constraints (in GramSchmidt)\n" 348 432 for rel in group: 349 errmsg += FormatConstraint(constrDict[rel],fixedList[rel])433 errmsg += _FormatConstraint(constrDict[rel],fixedList[rel]) 350 434 errmsg += "\n" 351 #import traceback352 #print traceback.format_exc()353 435 continue 354 355 436 mapvar = [] 356 437 group = group[:] … … 386 467 warnmsg += '\nPlease report this unexpected error\n' 387 468 for rel in group: 388 warnmsg += FormatConstraint(constrDict[rel],fixedList[rel])469 warnmsg += _FormatConstraint(constrDict[rel],fixedList[rel]) 389 470 warnmsg += "\n" 390 471 warnmsg += '\n\tNot refined: ' + notvaried + '\n' 391 472 try: 392 np.linalg.inv( constrArr)473 np.linalg.inv(multarr) 393 474 except: 394 475 errmsg += "\nSingular input. " … … 398 479 errmsg += 'Please report this unexpected error\n' 399 480 for rel in group: 400 errmsg += FormatConstraint(constrDict[rel],fixedList[rel])481 errmsg += _FormatConstraint(constrDict[rel],fixedList[rel]) 401 482 errmsg += "\n" 402 483 return errmsg,warnmsg … … 407 488 and stores them in global variables Also checks for internal 408 489 conflicts or inconsistencies in parameter/variable definitions. 409 Input: 410 groups,parmlist: see GroupConstraints 411 varyList: a list of parameters that will be varied 412 constrDict: a list of dicts defining relationships/constraints 413 constrDict = [{<constr1>}, {<constr2>}, ...] 414 where {<constr1>} is {'param1': mult1, 'param2': mult2,...} 415 fixedList: a list of values for each constraint/variable in 416 constrDict, value is either a float (contraint) or None (for 417 a new variable). 490 491 :param list groups: a list of grouped contraints where each constraint grouped containts a list of 492 indices for constraint constrDict entries, created in :func:`GroupConstraints` (returned as 1st value) 493 494 :param list parmlist: a list containing lists of parameter names contained in each group, 495 created in :func:`GroupConstraints` (returned as 1st value) 496 497 :param list varyList: a list of parameters names (strings of form ``<ph>:<hst>:<nam>``) that will be varied 498 499 :param dict constrDict: a list of dicts defining relationships/constraints (as defined in :func:`GroupConstraints`) 500 501 :param list fixedList: a list of values specifying a fixed value for each dict in constrDict. Values are 502 either strings that can be converted to floats, float values or None if the constraint defines a new 503 parameter 504 505 :param dict constrDict: a list of dicts defining relationships/constraints 506 418 507 ''' 419 508 global dependentParmList,arrayList,invarrayList,indParmList,consNum … … 424 513 if len(cdict) == 1: 425 514 fixedVarList.append(cdict.keys()[0]) 426 #print 'fixedVarList',fixedVarList427 515 428 516 # process equivalences: make a list of dependent and independent vars … … 525 613 if var in fixedVarList: 526 614 msg += '\nError: parameter '+var+" is Fixed and used in a constraint:\n\t" 527 msg += FormatConstraint(constrDict[rel],fixedList[rel])+"\n"615 msg += _FormatConstraint(constrDict[rel],fixedList[rel])+"\n" 528 616 # if var in equivVarList: 529 617 # msg += '\nError: parameter '+var+" is Equivalenced and used in a constraint:\n\t" 530 # msg += FormatConstraint(constrDict[rel],fixedList[rel])+"\n"618 # msg += _FormatConstraint(constrDict[rel],fixedList[rel])+"\n" 531 619 if varied > 0 and varied != len(constrDict[rel]): 532 620 msg += "\nNot all variables refined in constraint:\n\t" 533 msg += FormatConstraint(constrDict[rel],fixedList[rel])621 msg += _FormatConstraint(constrDict[rel],fixedList[rel]) 534 622 msg += '\nNot refined: ' + notvaried + '\n' 535 623 if fixedList[rel] is not None and varied > 0: … … 545 633 for group,varlist in zip(groups,parmlist): 546 634 if len(varlist) == 1: continue 547 arr = MakeArray(constrDict, group, varlist) 548 constrArr = FillInMissingRelations(arr,len(group)) 635 if len(varlist) < len(group): # too many relationships  no can do 636 msg = 'too many relationships' 637 try: 638 arr = _FillArray(group,constrDict,varlist) 639 _RowEchelon(len(group),arr,varlist) 640 constrArr = _FillArray(group,constrDict,varlist,FillDiagonals=True) 641 GramSchmidtOrtho(constrArr,len(group)) 642 except: 643 msg = 'Singular relationships' 644 549 645 mapvar = [] 550 646 group = group[:] … … 613 709 '''Takes a list of dependent parameter(s) and stores their 614 710 relationship to a single independent parameter (independentVar) 615 input: 616 independentVar  name of parameter that will set others (may 617 be varied) 618 dependentList  list of parameter that will set from 619 independentVar (may not be varied). Each item can be a parameter 711 712 :param str independentVar: name of master parameter that will be used to determine the value 713 to set the dependent variables 714 715 :param list dependentList: a list of parameters that will set from 716 independentVar. Each item in the list can be a string with the parameter 620 717 name or a tuple containing a name and multiplier: 621 ('parm1',('parm2',.5),) 718 ``['parm1',('parm2',.5),]`` 719 622 720 ''' 623 721 … … 643 741 def GetDependentVars(): 644 742 '''Return a list of dependent variables: e.g. variables that are 645 constrained in terms of other variables''' 743 constrained in terms of other variables 744 745 :returns: a list of variable names 746 747 ''' 646 748 dependentVars = [] 647 749 global dependentParmList … … 652 754 def GetIndependentVars(): 653 755 '''Return a list of independent variables: e.g. variables that are 654 created by constraints of other variables''' 756 created by constraints of other variables 757 758 :returns: a list of variable names 759 760 ''' 655 761 independentVars = [] 656 762 global indParmList,fixedDict … … 717 823 return sigmaDict 718 824 719 def FormatConstraint(RelDict,RelVal):825 def _FormatConstraint(RelDict,RelVal): 720 826 '''Formats a Constraint or Function for use in a convenient way''' 721 827 linelen = 45 … … 742 848 743 849 def VarRemapShow(varyList): 744 '''List out the saved relationships. 745 Returns a string containing the details. 850 '''List out the saved relationships. This should be done after the constraints have been 851 defined using :func:`StoreEquivalence`, :func:`GroupConstraints` and :func:`GenerateConstraints`. 852 853 :returns: a string containing the details of the contraint relationships 746 854 ''' 747 855 s = '' … … 792 900 '''Compute derivatives for Independent Parameters from the 793 901 derivatives for the original parameters 902 903 :param list varyList: a list of parameters names that will be varied 904 905 :param dict derivDict: a dict containing derivatives for parameter values keyed by the 906 parameter names. 907 908 :param dict dMdv: a dict containing derivatives for dependent parameter computed from 909 derivDict 910 794 911 ''' 795 912 global dependentParmList,arrayList,invarrayList,indParmList,invarrayList … … 816 933 Removes dependent variables from the varyList 817 934 818 This should be done once, before any variable refinement is done 819 to complete the parameter dictionary with the Independent Parameters 935 This should be done once, after the constraints have been 936 defined using :func:`StoreEquivalence`, 937 :func:`GroupConstraints` and :func:`GenerateConstraints` and 938 before any variable refinement is done 939 to complete the parameter dictionary by defining independent 940 parameters and satisfying the constraint equations. 941 942 :param dict parmDict: a dict containing parameter values keyed by the 943 parameter names. 944 This will contain updated values for both dependent and independent 945 parameters after Dict2Map is called. It will also contain some 946 unexpected entries of every constant value {'0':0.0} & {'1.0':1.0}, 947 which do not cause any problems. 948 949 :param list varyList: a list of parameters names that will be varied 950 951 820 952 ''' 821 953 # process the Independent vars: remove dependent ones from varylist … … 844 976 845 977 def Dict2Map(parmDict,varyList): 846 '''Convert the remapped values back to the original parameters 978 '''Applies the constraints defined using :func:`StoreEquivalence`, 979 :func:`GroupConstraints` and :func:`GenerateConstraints` by changing 980 values in a dict containing the parameters. This should be 981 done before the parameters are used for any computations 982 983 :param dict parmDict: a dict containing parameter values keyed by the 984 parameter names. 985 This will contain updated values for both dependent and independent 986 parameters after Dict2Map is called. It will also contain some 987 unexpected entries of every constant value {'0':0.0} & {'1.0':1.0}, 988 which do not cause any problems. 989 990 :param list varyList: a list of parameters names that will be varied 847 991 848 This should be done to apply constraints to parameter values (use 849 Map2Dict once first). It should be done as part of the Model function 850 evaluation, before any computation is done 851 ''' 852 #I think this needs fixing to update parmDict with new values 853 # from the constraints based on what is in varyList  RVD. Don't think so  BHT 992 ''' 854 993 global dependentParmList,arrayList,invarrayList,indParmList,fixedDict 855 994 # reset fixed values (should not be needed, but very quick) … … 867 1006 # internal routines follow (these routines are unlikely to be called 868 1007 # from outside the module) 869 def FillInMissingRelations(arrayin,nvars): 870 '''Fill in unspecified rows in array with noncolinear unit vectors 871 arrayin is a square array, where only the first nvars rows are defined. 872 873 Returns a new array where all rows are defined 874 875 Throws an exception if there is no nonsignular result 876 (likely because two or more input rows are colinear) 877 ''' 878 a = arrayin.copy() 879 n = nvars 880 # fill in the matrix with basis vectors not colinear with any of the starting set 881 for i in range(n,len(a)): 882 try: 883 a[n] = VectorProduct(a[:n]) 884 except: 885 raise Exception,"VectorProduct failed, singular input?" 886 n += 1 887 # use the GS algorithm to compute a complete set of unit vectors orthogonal 888 # to the first in array 889 gs = GramSchmidtOrtho(a) 890 n = nvars 891 # now replace the created vectors with the ones least correlated to the 892 # initial set 893 vlist = [v for v in gs[1:]] # drop the first row 894 for j in range(n,len(a)): 895 minavg = None 896 droplist = [] 897 for k in range(len(vlist)): 898 v = vlist[k] 899 avgcor = 0 900 for i in range(n): 901 cor = np.dot(a[i],vlist[k])/(np.linalg.norm(a[i]) * np.linalg.norm(vlist[k]) ) 902 if np.allclose(cor, 1.0): 903 droplist.append(k) # found a vector colinear to one we already have 904 # don't need it again 905 #vlist.pop(k) 906 break 907 avgcor += cor 908 else: 909 if minavg == None: 910 minavg = abs(avgcor/n) 911 minnum = k 912 elif abs(avgcor/n) < minavg: 913 minavg = abs(avgcor/n) 914 minnum = k 915 if minavg == None: 916 raise Exception,"Failed to find a noncolinear GS vector for row %d. Should not happen!" % n 917 a[j] = vlist[minnum] 918 droplist.append(minnum) # don't need to use this vector again 919 for i in sorted(droplist,reverse=True): 920 vlist.pop(i) 921 n += 1 922 return a 923 924 925 def MakeArray(constrDict, group, varlist): 926 """Convert the multipliers in a constraint group to an array of 927 coefficients and place in a square numpy array, adding empty rows as 928 needed. 929 930 Throws an exception if some sanity checks fail. 931 """ 932 # do some error checks 933 if len(varlist) < len(group): # too many relationships  no can do 934 raise Exception,'The number of relationships (%d) exceeds the number of parameters (%d):\n\t%s\n\t%s'% ( 935 len(group),len(varlist),group,varlist) 936 # put the coefficients into an array 937 multarr = np.zeros(2*[len(varlist),]) 938 i = 0 939 for cnum in group: 940 cdict = constrDict[cnum] 941 j = 0 942 for var in varlist: 943 m = cdict.get(var) 944 if m != None: 945 multarr[i,j] = m 946 j += 1 947 i += 1 948 return multarr 949 950 def GramSchmidtOrtho(arrayin): 1008 1009 def GramSchmidtOrtho(a,nkeep=0): 951 1010 '''Use the GramSchmidt process (http://en.wikipedia.org/wiki/GramSchmidt) to 952 1011 find orthonormal unit vectors relative to first row. 1012 1013 If nkeep is nonzero, the first nkeep rows in the array are not changed 1014 953 1015 input: 954 arrayin: a 2D nonsi gnular square array1016 arrayin: a 2D nonsingular square array 955 1017 returns: 956 1018 a orthonormal set of unit vectors as a square array … … 959 1021 'Projection operator' 960 1022 return a*(np.dot(a,b)/np.dot(a,a)) 961 a = arrayin.copy() 962 for j in range (len(a)): 1023 for j in range(nkeep,len(a)): 963 1024 for i in range(j): 964 a[j] = a[j] proj(a[i],a[j])1025 a[j] = proj(a[i],a[j]) 965 1026 if np.allclose(np.linalg.norm(a[j]),0.0): 966 1027 raise Exception,"Singular input to GramSchmidtOrtho" 967 a[j] = a[j]/np.linalg.norm(a[j])1028 a[j] /= np.linalg.norm(a[j]) 968 1029 return a 969 1030 970 def VectorProduct(vl): 971 '''Evaluate the ndimensional "cross" product of the list of vectors in vl 972 vl can be any length. The length of each vector is arbitrary, but 973 all must be the same length. See http://en.wikipedia.org/wiki/LeviCivita_symbol 974 975 This appears to return a vector perpendicular to the supplied set. 976 977 Throws an exception if a vector can not be obtained because the input set of 978 vectors is already complete or contains any redundancies. 979 980 Uses LeviCitvitaMult recursively to obtain all permutations of vector elements 981 ''' 982 i = 0 983 newvec = [] 984 for e in vl[0]: 985 i += 1 986 tl = [([i,],1),] 987 seq = LeviCitvitaMult(tl ,vl) 988 tsum = 0 989 for terms,prod in seq: 990 tsum += EvalLCterm(terms) * prod 991 newvec.append(tsum) 992 if np.allclose(newvec,0.0): 993 raise Exception,"VectorProduct failed, singular or already complete input" 994 return newvec 995 996 def LeviCitvitaMult(tin,vl): 997 '''Recursion formula to compute all permutations of vector element products 998 The first term in each list is the indices of the LeviCivita symbol, note 999 that if an index is repeated, the value is zero, so the element is not included. 1000 The second term is the product of the vector elements 1001 ''' 1002 v = vl[0] 1003 vl1 = vl[1:] 1004 1005 j = 0 1006 tl = [] 1007 for e in vl[0]: 1008 j += 1 1009 for ind,vals in tin: 1010 if j in ind: continue 1011 tl.append((ind+[j,],vals*e)) 1012 if len(vl1): 1013 return LeviCitvitaMult(tl,vl1) 1031 def _FillArray(sel,dict,collist,FillDiagonals=False): 1032 '''Construct a n by n matrix (n = len(collist) 1033 filling in the rows using the relationships defined 1034 in the dictionaries selected by sel 1035 1036 If FillDiagonals is True, diagonal elements in the 1037 array are set to 1.0 1038 ''' 1039 n = len(collist) 1040 if FillDiagonals: 1041 arr = np.eye(n) 1014 1042 else: 1015 return tl 1016 1017 def EvalLCterm(term): 1018 '''Evaluate the LeviCivita symbol Epsilon(i,j,k,...)''' 1019 p = 1 1020 for i in range(len(term)): 1021 for j in range(i+1,len(term)): 1022 p *= (term[j]  term[i]) 1023 if p == 0: return 0 1024 return p/abs(p) 1025 1043 arr = np.zeros(2*[n,]) 1044 # fill the top rows 1045 for i,cnum in enumerate(sel): 1046 for j,var in enumerate(collist): 1047 arr[i,j] = dict[cnum].get(var,0) 1048 return arr 1049 1050 def _SwapColumns(i,m,v): 1051 '''Swap columns in matrix m as well as the labels in v 1052 so that element (i,i) is replaced by the first nonzero element in row i after that element 1053 1054 Throws an exception if there are no nonzero elements in that row 1055 ''' 1056 for j in range(i+1,len(v)): 1057 if not np.allclose(m[i,j],0): 1058 m[:,(i,j)] = m[:,(j,i)] 1059 v[i],v[j] = v[j],v[i] 1060 return 1061 else: 1062 raise Exception,'Singular input' 1063 1064 def _RowEchelon(m,arr,collist): 1065 '''Convert the first m rows in Matrix arr to rowechelon form 1066 exchanging columns in the matrix and collist as needed. 1067 1068 throws an exception if the matrix is singular because 1069 the first m rows are not linearly independent 1070 ''' 1071 n = len(collist) 1072 for i in range(m): 1073 if np.allclose(arr[i,i],0): 1074 _SwapColumns(i,arr,collist) 1075 arr[i,:] /= arr[i,i] # normalize row 1076 # subtract current row from subsequent rows to set values to left of diagonal to 0 1077 for j in range(i+1,m): 1078 arr[j,:] = arr[i,:] * arr[j,i] 1026 1079 1027 1080 if __name__ == "__main__": 1028 import sys1029 1081 parmdict = {} 1030 1082 constrDict = [ … … 1047 1099 '2::atomy:3', '2::atomz:3', 1048 1100 '0:12:Scale', '0:11:Scale', '0:14:Scale', '0:13:Scale', '0:0:Scale'] 1049 varyList = ['0::A0', '0::AUiso:0', '0::Afrac:1', '0::Afrac:2', '0::Afrac:3', '0::Afrac:4', '0::dAx:5', '0::dAy:5', '0::dAz:5', '0::AUiso:5', ':0:Back:0', ':0:Back:1', ':0:Back:2', ':0:Back:3', ':0:Back:4', ':0:Back:5', ':0:Back:6', ':0:Back:7', ':0:Back:8', ':0:Back:9', ':0:Back:10', ':0:Back:11', ':0:U', ':0:V', ':0:W', ':0:X', ':0:Y', ':0:Scale', ':0:DisplaceX', ':0:DisplaceY'] 1050 constrDict = [ 1051 {'0::Afrac:4': 24.0, '0::Afrac:1': 16.0, '0::Afrac:3': 24.0, '0::Afrac:2': 16.0}, 1052 {'0::Afrac:1': 1.0, '0::Afrac:2': 1.0}, 1053 {'0::Afrac:4': 1.0, '0::Afrac:3': 1.0}] 1054 fixedList = ['40.0', '1.0', '1.0'] 1101 # e,w = CheckConstraints([, 1102 # [{'2:0:Scale': 1.0, '5:0:Scale': 1.0, '10:0:Scale': 1.0, '6:0:Scale': 1.0, '9:0:Scale': 1.0, '8:0:Scale': 1.0,# '3:0:Scale': 1.0, '4:0:Scale': 1.0, '7:0:Scale': 1.0, '1:0:Scale': 1.0, '0:0:Scale': 1.0}], 1103 # ['1.0']) 1104 # if e: print 'error=',e 1105 # if w: print 'error=',w 1106 # varyList = ['0::A0', '0::AUiso:0', '0::Afrac:1', '0::Afrac:2', '0::Afrac:3', '0::Afrac:4', '0::dAx:5', '0::dAy:5', '0::dAz:5', '0::AUiso:5', ':0:Back:0', ':0:Back:1', ':0:Back:2', ':0:Back:3', ':0:Back:4', ':0:Back:5', ':0:Back:6', ':0:Back:7', ':0:Back:8', ':0:Back:9', ':0:Back:10', ':0:Back:11', ':0:U', ':0:V', ':0:W', ':0:X', ':0:Y', ':0:Scale', ':0:DisplaceX', ':0:DisplaceY'] 1107 # constrDict = [ 1108 # {'0::Afrac:4': 24.0, '0::Afrac:1': 16.0, '0::Afrac:3': 24.0, '0::Afrac:2': 16.0}, 1109 # {'0::Afrac:1': 1.0, '0::Afrac:2': 1.0}, 1110 # {'0::Afrac:4': 1.0, '0::Afrac:3': 1.0}] 1111 # fixedList = ['40.0', '1.0', '1.0'] 1055 1112 1056 1113 errmsg, warnmsg = CheckConstraints(varylist,constrDict,fixedList) … … 1090 1147 for key in sorted(parmdict.keys()): 1091 1148 print ' '+key,'\t',before.get(key),'\t',parmdict[key] 1092 dMdv = len(varylist)*[0]1093 deriv = {}1094 for i,v in enumerate(parmdict.keys()): deriv[v]=i1095 Dict2Deriv(varylist,deriv,dMdv)1149 # dMdv = len(varylist)*[0] 1150 # deriv = {} 1151 # for i,v in enumerate(parmdict.keys()): deriv[v]=i 1152 # Dict2Deriv(varylist,deriv,dMdv)
Note: See TracChangeset
for help on using the changeset viewer.