Changeset 3041


Ignore:
Timestamp:
Sep 3, 2017 11:18:33 AM (6 years ago)
Author:
toby
Message:

Fully implement multiprocessing on reflection loops, but currently disabled, when enabled default is still not to use unless turned on in config; add timing code

Location:
trunk
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIdataGUI.py

    r3037 r3041  
    27652765        arg = sys.argv
    27662766        if len(arg) > 1 and arg[1]:
    2767             self.GSASprojectfile = os.path.splitext(arg[1])[0]+'.gpx'
     2767            try:
     2768                self.GSASprojectfile = os.path.splitext(arg[1])[0]+u'.gpx'
     2769            except:
     2770                self.GSASprojectfile = os.path.splitext(arg[1])[0]+'.gpx'
    27682771            self.dirname = os.path.abspath(os.path.dirname(arg[1]))
    27692772            if self.dirname: os.chdir(self.dirname)
    … …  
    36073610            GetGPX()
    36083611        else:
    3609             self.GSASprojectfile = os.path.splitext(filename)[0]+'.gpx'
     3612            try:
     3613                self.GSASprojectfile = os.path.splitext(filename)[0]+u'.gpx'
     3614            except:
     3615                self.GSASprojectfile = os.path.splitext(filename)[0]+'.gpx'
    36103616            self.dirname = os.path.split(filename)[0]
    36113617

  • trunk/GSASIImpsubs.py

    r3000 r3041  
    4444
    4545def InitMP(allowMP=True):
    46     '''Called in routines to initialize use of Multiprocessing
     46    '''Called to initialize use of Multiprocessing
    4747    '''
    4848    global useMP,ncores
    49     if ncores is not None: return
     49    if ncores is not None: return useMP,ncores
    5050    useMP = False
    5151    if not allowMP:
    5252        print('Multiprocessing disabled')
    5353        ncores = 0
    54         return
    55     ncores = GSASIIpath.GetConfigValue('Multiprocessing_cores',-1)
    56     if ncores < 0: ncores = mp.cpu_count()
     54        return useMP,ncores
     55    ncores = GSASIIpath.GetConfigValue('Multiprocessing_cores',0)
     56    if ncores < 0: ncores = mp.cpu_count()/2
    5757    if ncores > 1:
    5858        useMP = True
    59     #if GSASIIpath.GetConfigValue('debug'):
    60     if True:
     59    #if GSASIIpath.GetConfigValue('debug') and useMP:
     60    if useMP:
    6161        print('Multiprocessing with {} cores enabled'.format(ncores))
    62 
    63 ################################################################################
    64 # derivative computation
    65 ################################################################################       
    66 def InitDerivGlobals(im1,calcControls1,SGMT1,hfx1,phfx1,pfx1,G1,GB1,g1,SGData1,
    67                      parmDict1,wave1,shl1,x1,cw1,Ka21,A1,varylist1,dependentVars1,
    68                      dFdvDict1,lamRatio1,kRatio1,doPawley1,pawleyLookup1):
    69     '''Initialize for the computation of derivatives. Puts lots of junk into the global
    70     namespace in this module, including the arrays for derivatives (when needed.)
    71     '''
    72     global im,calcControls,SGMT,hfx,phfx,pfx,G,GB,g,SGData,parmDict,wave,shl,x,cw,Ka2,A
    73     global varylist,dependentVars,dFdvDict,lamRatio,kRatio,doPawley,pawleyLookup
    74     im = im1
    75     calcControls = calcControls1
    76     SGMT = SGMT1
    77     hfx = hfx1
    78     phfx = phfx1
    79     pfx = pfx1
    80     G = G1
    81     GB = GB1
    82     g = g1
    83     SGData = SGData1
    84     parmDict = parmDict1
    85     wave = wave1
    86     shl = shl1
    87     x = ma.getdata(x1)
    88     cw = cw1
    89     Ka2 = Ka21
    90     A = A1
    91     varylist = varylist1
    92     dependentVars = dependentVars1
    93     dFdvDict = dFdvDict1
    94     lamRatio = lamRatio1
    95     kRatio = kRatio1
    96     doPawley = doPawley1
    97     pawleyLookup = pawleyLookup1
    98     # determine the parameters that will have derivatives computed only at end
    99     global nonatomvarylist
    100     nonatomvarylist = []
    101     for name in varylist:
    102         if '::RBV;' not in name:
    103             try:
    104                 aname = name.split(pfx)[1][:2]
    105                 if aname not in ['Af','dA','AU','RB','AM','Xs','Xc','Ys','Yc','Zs','Zc',    \
    106                     'Tm','Xm','Ym','Zm','U1','U2','U3']: continue # skip anything not an atom or rigid body param
    107             except IndexError:
    108                 continue
    109         nonatomvarylist.append(name)
    110     global nonatomdependentVars
    111     nonatomdependentVars = []
    112     for name in dependentVars:
    113         if '::RBV;' not in name:
    114             try:
    115                 aname = name.split(pfx)[1][:2]
    116                 if aname not in ['Af','dA','AU','RB','AM','Xs','Xc','Ys','Yc','Zs','Zc',    \
    117                     'Tm','Xm','Ym','Zm','U1','U2','U3']: continue # skip anything not an atom or rigid body param
    118             except IndexError:
    119                 continue
    120         nonatomdependentVars.append(name)
    121     # create local copies of derivative arrays, if multiprocessing will be used
    122     if useMP:
    123         global dMdv
    124         dMdv = np.zeros(shape=(len(varylist),len(x)))
    125         global depDerivDict
    126         depDerivDict = {j:np.zeros(shape=len(x)) for j in dependentVars}
    127            
    128 
    129 def cellVaryDerv(pfx,SGData,dpdA):
    130     if SGData['SGLaue'] in ['-1',]:
    131         return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
    132             [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
    133     elif SGData['SGLaue'] in ['2/m',]:
    134         if SGData['SGUniq'] == 'a':
    135             return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
    136         elif SGData['SGUniq'] == 'b':
    137             return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
    138         else:
    139             return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
    140     elif SGData['SGLaue'] in ['mmm',]:
    141         return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
    142     elif SGData['SGLaue'] in ['4/m','4/mmm']:
    143         return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
    144     elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
    145         return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
    146     elif SGData['SGLaue'] in ['3R', '3mR']:
    147         return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
    148     elif SGData['SGLaue'] in ['m3m','m3']:
    149         return [[pfx+'A0',dpdA[0]]]
    150 
    151 def ComputeDerivMPbatch(reflsList):
    152     '''Computes a the derivatives for a batch of reflections and sums the them into
    153     global arrays dMdv & depDerivDict. These arrays are returned once the computation
    154     is completed.
    155     '''
    156     for refl,iref,fmin,fmax,iBeg,iFin in reflsList:
    157         if ComputeDeriv(refl,iref,fmin,fmax,iBeg,iFin,dMdv,depDerivDict): break
    158     return dMdv,depDerivDict
    159 
    160 def ComputeDeriv(refl,iref,fmin,fmax,iBeg,iFin,dMdv,depDerivDict):
    161     '''Compute the parameter derivatives for a single reflection and add the results
    162     into either array dMdv or depDerivDict
    163     '''
    164     global wave
    165     if im:
    166         h,k,l,m = refl[:4]
    167     else:
    168         h,k,l = refl[:3]
    169     Uniq = np.inner(refl[:3],SGMT)
    170     if 'T' in calcControls[hfx+'histType']:
    171         wave = refl[14+im]
    172     dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = G2stMth.GetIntensityDerv(refl,im,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
    173     pos = refl[5+im]
    174     calcKa2 = False
    175     if 'C' in calcControls[hfx+'histType']:
    176         tanth = tand(pos/2.0)
    177         costh = cosd(pos/2.0)
    178         if Ka2:
    179             pos2 = refl[5+im]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
    180             iBeg2 = np.searchsorted(x,pos2-fmin)
    181             iFin2 = np.searchsorted(x,pos2+fmax)
    182             if iBeg2-iFin2:
    183                 calcKa2 = True
    184                 iFin = iFin2       
    185         lenBF = iFin-iBeg
    186         dMdpk = np.zeros(shape=(6,lenBF))
    187         dMdipk = G2pwd.getdFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,x[iBeg:iFin])
    188         for i in range(5):
    189             dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11+im]*refl[9+im]*dMdipk[i]
    190         dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
    191         if calcKa2:
    192             dMdpk2 = np.zeros(shape=(6,lenBF))
    193             dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,x[iBeg:iFin])
    194             for i in range(5):
    195                 dMdpk2[i] = 100.*cw[iBeg:iFin]*refl[11+im]*refl[9+im]*kRatio*dMdipk2[i]
    196             dMdpk2[5] = 100.*cw[iBeg:iFin]*refl[11+im]*dMdipk2[0]
    197             dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
    198     else:   #'T'OF
    199         lenBF = iFin-iBeg
    200         if lenBF < 0: return True  #bad peak coeff
    201         dMdpk = np.zeros(shape=(6,lenBF))
    202         dMdipk = G2pwd.getdEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],x[iBeg:iFin])
    203         for i in range(6):
    204             dMdpk[i] += refl[11+im]*refl[9+im]*dMdipk[i]      #cw[iBeg:iFin]*
    205         dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}           
    206     if doPawley:
    207         try:
    208             if im:
    209                 pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
    210             else:
    211                 pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
    212             idx = varylist.index(pIdx)
    213             dMdv[idx][iBeg:iFin] = dervDict['int']/refl[9+im]
    214             if Ka2: #not for TOF either
    215                 dMdv[idx][iBeg:iFin] += dervDict2['int']/refl[9+im]
    216         except: # ValueError:
    217             pass
    218     if 'C' in calcControls[hfx+'histType']:
    219         dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY,dpdV = G2stMth.GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict)
    220         names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
    221             hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
    222             hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
    223             hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
    224             hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
    225             hfx+'DisplaceY':[dpdY,'pos'],}
    226         if 'Bragg' in calcControls[hfx+'instType']:
    227             names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
    228                 hfx+'SurfRoughB':[dFdAb[1],'int'],})
    229         else:
    230             names.update({hfx+'Absorption':[dFdAb,'int'],})
    231     else:   #'T'OF
    232         dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV = G2stMth.GetReflPosDerv(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)
    233         names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
    234             hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
    235             hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4+im],'gam'],hfx+'Y':[refl[4+im]**2,'gam'],
    236             hfx+'alpha':[1./refl[4+im],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4+im]**4,'bet'],
    237             hfx+'beta-q':[1./refl[4+im]**2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4+im]**2,'sig'],
    238             hfx+'sig-2':[refl[4+im]**4,'sig'],hfx+'sig-q':[1./refl[4+im]**2,'sig'],
    239             hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
    240     for name in names:
    241         item = names[name]
    242         if name in varylist:
    243             dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
    244             if calcKa2:
    245                 dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict2[item[1]]
    246         elif name in dependentVars:
    247             depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
    248             if calcKa2:
    249                 depDerivDict[name][iBeg:iFin] += item[0]*dervDict2[item[1]]
    250     for iPO in dIdPO:
    251         if iPO in varylist:
    252             dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
    253             if calcKa2:
    254                 dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict2['int']
    255         elif iPO in dependentVars:
    256             depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
    257             if calcKa2:
    258                 depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict2['int']
    259     for i,name in enumerate(['omega','chi','phi']):
    260         aname = pfx+'SH '+name
    261         if aname in varylist:
    262             dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
    263             if calcKa2:
    264                 dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict2['int']
    265         elif aname in dependentVars:
    266             depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
    267             if calcKa2:
    268                 depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict2['int']
    269     for iSH in dFdODF:
    270         if iSH in varylist:
    271             dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
    272             if calcKa2:
    273                 dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict2['int']
    274         elif iSH in dependentVars:
    275             depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
    276             if calcKa2:
    277                 depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict2['int']
    278     cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
    279     for name,dpdA in cellDervNames:
    280         if name in varylist:
    281             dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
    282             if calcKa2:
    283                 dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict2['pos']
    284         elif name in dependentVars: #need to scale for mixed phase constraints?
    285             depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
    286             if calcKa2:
    287                 depDerivDict[name][iBeg:iFin] += dpdA*dervDict2['pos']
    288     dDijDict = G2stMth.GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict)
    289     for name in dDijDict:
    290         if name in varylist:
    291             dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
    292             if calcKa2:
    293                 dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict2['pos']
    294         elif name in dependentVars:
    295             depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
    296             if calcKa2:
    297                 depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict2['pos']
    298     for i,name in enumerate([pfx+'mV0',pfx+'mV1',pfx+'mV2']):
    299         if name in varylist:
    300             dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict['pos']
    301             if calcKa2:
    302                 dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict2['pos']
    303         elif name in dependentVars:
    304             depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict['pos']
    305             if calcKa2:
    306                 depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict2['pos']
    307     if 'C' in calcControls[hfx+'histType']:
    308         sigDict,gamDict = G2stMth.GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
    309     else:   #'T'OF
    310         sigDict,gamDict = G2stMth.GetSampleSigGamDerv(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
    311     for name in gamDict:
    312         if name in varylist:
    313             dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
    314             if calcKa2:
    315                 dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict2['gam']
    316         elif name in dependentVars:
    317             depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
    318             if calcKa2:
    319                 depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict2['gam']
    320     for name in sigDict:
    321         if name in varylist:
    322             dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
    323             if calcKa2:
    324                 dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict2['sig']
    325         elif name in dependentVars:
    326             depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
    327             if calcKa2:
    328                 depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict2['sig']
    329     for name in ['BabA','BabU']:
    330         if refl[9+im]:
    331             if phfx+name in varylist:
    332                 dMdv[varylist.index(phfx+name)][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict['int']/refl[9+im]
    333                 if calcKa2:
    334                     dMdv[varylist.index(phfx+name)][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict2['int']/refl[9+im]
    335             elif phfx+name in dependentVars:                   
    336                 depDerivDict[phfx+name][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict['int']/refl[9+im]
    337                 if calcKa2:
    338                     depDerivDict[phfx+name][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict2['int']/refl[9+im]                 
    339     if not doPawley and not parmDict[phfx+'LeBail']:
    340         #do atom derivatives -  for RB,F,X & U so far - how do I scale mixed phase constraints?
    341         corr = 0.
    342         #corr2 = 0.
    343         if refl[9+im]:             
    344             corr = dervDict['int']/refl[9+im]
    345             #if calcKa2:  # commented out in Bob's code. Why?
    346             #    corr2 = dervDict2['int']/refl[9+im]
    347         for name in nonatomvarylist:
    348             dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
    349             #if calcKa2:
    350             #   dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr2 # unneeded w/o above
    351         for name in nonatomdependentVars:
    352            depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
    353            #if calcKa2:
    354            #    depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr2
     62    return useMP,ncores
    35563
    35664################################################################################
    35765# Fobs Squared computation
    35866################################################################################       
    359 #x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2
    36067def InitFobsSqGlobals(x1,ratio1,shl1,xB1,xF1,im1,lamRatio1,kRatio1,xMask1,Ka21):
    36168    '''Initialize for the computation of Fobs Squared for powder histograms.
    … …  
    464171    cw = np.append(cw,cw[-1])
    465172    # create local copies of ycalc array
    466     if useMP:
    467         global yc
    468         yc = np.zeros_like(x1)
     173    global yc
     174    yc = np.zeros_like(x)
    469175
    470176

  • trunk/GSASIIobj.py

    r3023 r3041  
    27072707    return Id
    27082708
     2709class ShowTiming(object):
     2710    '''An object to use for timing repeated sections of code.
     2711
     2712    Create the object with::
     2713       tim0 = ShowTiming()
     2714
     2715    Tag sections of code to be timed with::
     2716       tim0.start('start')
     2717       tim0.start('in section 1')
     2718       tim0.start('in section 2')
     2719    etc. (Note that each section should have a unique label.)
     2720
     2721    After the last section, end timing with::
     2722       tim0.end()
     2723
     2724    Show timing results with::
     2725       tim0.show()
     2726       
     2727    '''
     2728    def __init__(self):
     2729        self.timeSum =  []
     2730        self.timeStart = []
     2731        self.label = []
     2732        self.prev = None
     2733    def start(self,label):
     2734        if label in self.label:
     2735            i = self.label.index(label)
     2736            self.timeStart[i] = time.time()
     2737        else:
     2738            i = len(self.label)
     2739            self.timeSum.append(0.0)
     2740            self.timeStart.append(time.time())
     2741            self.label.append(label)
     2742        if self.prev is not None:
     2743            self.timeSum[self.prev] += self.timeStart[i] - self.timeStart[self.prev]
     2744        self.prev = i
     2745    def end(self):
     2746        if self.prev is not None:
     2747            self.timeSum[self.prev] += time.time() - self.timeStart[self.prev]
     2748        self.prev = None
     2749    def show(self):
     2750        sumT = sum(self.timeSum)
     2751        print('Timing results (total={:.2f} sec)'.format(sumT))
     2752        for i,(lbl,val) in enumerate(zip(self.label,self.timeSum)):
     2753            print('{} {:20} {:8.2f} ms {:5.2f}%'.format(i,lbl,1000.*val,100*val/sumT))
     2754
    27092755
    27102756if __name__ == "__main__":

  • trunk/GSASIIstrMain.py

    r3023 r3041  
    771771
    772772if __name__ == '__main__':
     773    GSASIIpath.InvokeDebugOpts()
    773774    main()

  • trunk/GSASIIstrMath.py

    r3005 r3041  
    1717import numpy.linalg as nl
    1818import scipy.stats as st
     19import multiprocessing as mp
    1920import GSASIIpath
    2021GSASIIpath.SetVersionNumber("$Revision$")
    … …  
    2627import GSASIImath as G2mth
    2728import GSASIIobj as G2obj
     29import GSASIImpsubs as G2mp
     30G2mp.InitMP(False)  # This disables multiprocessing
    2831
    2932sind = lambda x: np.sin(x*np.pi/180.)
    … …  
    29232926               
    29242927def GetFobsSq(Histograms,Phases,parmDict,calcControls):
    2925     'Compute the observed structure factors for Powder histograms'
    2926     #starttime = time.time(); print 'start GetFobsSq'
     2928    '''Compute the observed structure factors for Powder histograms and store in reflection array
     2929    Multiprocessing support added
     2930    '''
     2931    if GSASIIpath.GetConfigValue('debug'):
     2932        starttime = time.time() #; print 'start GetFobsSq'
    29272933    histoList = Histograms.keys()
    29282934    histoList.sort()
     2935    Ka2 = shl = lamRatio = kRatio = None
    29292936    for histogram in histoList:
    29302937        if 'PWDR' in histogram[:4]:
    … …  
    29662973                sumInt = 0.0
    29672974                nExcl = 0
    2968                 for refl in refDict['RefList']:
    2969                     if 'C' in calcControls[hfx+'histType']:
    2970                         yp = np.zeros_like(yb)
    2971                         Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
    2972                         iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
    2973                         iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
    2974                         iFin2 = iFin
    2975                         if not iBeg+iFin:       #peak below low limit - skip peak
    2976                             continue
    2977                         if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
    2978                             refl[3+im] *= -1
    2979                             nExcl += 1
    2980                             continue
    2981                         elif not iBeg-iFin:     #peak above high limit - done
    2982                             break
    2983                         elif iBeg < iFin:
    2984                             yp[iBeg:iFin] = refl[11+im]*refl[9+im]*G2pwd.getFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))    #>90% of time spent here
    2985                             sumInt += refl[11+im]*refl[9+im]
    2986                             if Ka2:
    2987                                 pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
    2988                                 Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6+im],refl[7+im],shl)
    2989                                 iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
    2990                                 iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
    2991                                 if iFin2 > iBeg2:
    2992                                     yp[iBeg2:iFin2] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2]))        #and here
    2993                                     sumInt += refl[11+im]*refl[9+im]*kRatio
    2994                             refl[8+im] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11+im]*(1.+kRatio)),0.0))
    2995                             if parmDict[phfx+'LeBail']:
    2996                                 refl[9+im] = refl[8+im]
    2997                                
    2998                     elif 'T' in calcControls[hfx+'histType']:
    2999                         yp = np.zeros_like(yb)
    3000                         Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
    3001                         iBeg = max(xB,np.searchsorted(x,refl[5+im]-fmin))
    3002                         iFin = max(xB,min(np.searchsorted(x,refl[5+im]+fmax),xF))
    3003                         if not iBeg+iFin:       #peak below low limit - skip peak
    3004                             continue
    3005                         if ma.all(xMask[iBeg:iFin]):    #peak entirely masked - skip peak
    3006                             refl[3+im] *= -1
    3007                             nExcl += 1
    3008                             continue
    3009                         elif not iBeg-iFin:     #peak above high limit - done
    3010                             break
    3011                         if iBeg < iFin:
    3012                             yp[iBeg:iFin] = refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))  #>90% of time spent here
    3013                             refl[8+im] = np.sum(np.where(ratio[iBeg:iFin]>0.,yp[iBeg:iFin]*ratio[iBeg:iFin]/refl[11+im],0.0))
    3014                             if parmDict[phfx+'LeBail']:
    3015                                 refl[9+im] = refl[8+im]
    3016                             sumInt += refl[11+im]*refl[9+im]
     2975                # test to see if we are using multiprocessing below
     2976                useMP,ncores = G2mp.InitMP()
     2977                if len(refDict['RefList']) < 100: useMP = False       
     2978                if useMP: # multiprocessing: create a set of initialized Python processes
     2979                    MPpool = mp.Pool(G2mp.ncores,G2mp.InitFobsSqGlobals,
     2980                                    [x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2])
     2981                    profArgs = [[] for i in range(G2mp.ncores)]
     2982                else:
     2983                    G2mp.InitFobsSqGlobals(x,ratio,shl,xB,xF,im,lamRatio,kRatio,xMask,Ka2)
     2984                if 'C' in calcControls[hfx+'histType']:
     2985                    # are we multiprocessing?
     2986                    for iref,refl in enumerate(refDict['RefList']):
     2987                        if useMP:
     2988                            profArgs[iref%G2mp.ncores].append((refl,iref))
     2989                        else:
     2990                            icod= G2mp.ComputeFobsSqCW(refl,iref)
     2991                            if type(icod) is tuple:
     2992                                refl[8+im] = icod[0]
     2993                                sumInt += icod[1]
     2994                                if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
     2995                            elif icod == -1:
     2996                                refl[3+im] *= -1
     2997                                nExcl += 1
     2998                            elif icod == -2:
     2999                                break
     3000                    if useMP:
     3001                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqCWbatch,profArgs):
     3002                            sumInt += sInt
     3003                            for refl8im,irefl in resList:
     3004                                if refl8im is None:
     3005                                    refDict['RefList'][irefl][3+im] *= -1
     3006                                    nExcl += 1
     3007                                else:
     3008                                    refDict['RefList'][irefl][8+im] = refl8im
     3009                                    if parmDict[phfx+'LeBail']:
     3010                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
     3011                        MPpool.terminate()
     3012                elif 'T' in calcControls[hfx+'histType']:
     3013                    for iref,refl in enumerate(refDict['RefList']):
     3014                        if useMP:
     3015                            profArgs[iref%G2mp.ncores].append((refl,iref))
     3016                        else:
     3017                            icod= G2mp.ComputeFobsSqTOF(refl,iref)
     3018                            if type(icod) is tuple:
     3019                                refl[8+im] = icod[0]
     3020                                sumInt += icod[1]
     3021                                if parmDict[phfx+'LeBail']: refl[9+im] = refl[8+im]
     3022                            elif icod == -1:
     3023                                refl[3+im] *= -1
     3024                                nExcl += 1
     3025                            elif icod == -2:
     3026                                break
     3027                    if useMP:
     3028                        for sInt,resList in MPpool.imap_unordered(G2mp.ComputeFobsSqTOFbatch,profArgs):
     3029                            sumInt += sInt
     3030                            for refl8im,irefl in resList:
     3031                                if refl8im is None:
     3032                                    refDict['RefList'][irefl][3+im] *= -1
     3033                                    nExcl += 1
     3034                                else:
     3035                                    refDict['RefList'][irefl][8+im] = refl8im
     3036                                    if parmDict[phfx+'LeBail']:
     3037                                        refDict['RefList'][irefl][9+im] = refDict['RefList'][irefl][8+im]
     3038                        MPpool.terminate()
     3039                sumFo = 0.0
     3040                sumdF = 0.0
     3041                sumFosq = 0.0
     3042                sumdFsq = 0.0
     3043                for iref,refl in enumerate(refDict['RefList']):
    30173044                    Fo = np.sqrt(np.abs(refl[8+im]))
    30183045                    Fc = np.sqrt(np.abs(refl[9]+im))
    … …  
    30333060            Histogram = Histograms[histogram]
    30343061            Histogram['Residuals']['hId'] = Histograms[histogram]['hId']
    3035     #print 'end GetFobsSq t=',time.time()-starttime
     3062    if GSASIIpath.GetConfigValue('debug'):
     3063        print 'GetFobsSq t=',time.time()-starttime
    30363064               
    30373065def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
    30383066    'Computes the powder pattern for a histogram based on contributions from all used phases'
    3039    
    3040     #starttime = time.time(); print 'start getPowderProfile'
     3067    if GSASIIpath.GetConfigValue('debug'): starttime = time.time()
    30413068   
    30423069    def GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict):
    … …  
    30863113        else:
    30873114            wave = parmDict[hfx+'Lam']
     3115    else:
     3116        shl = 0.
    30883117    for phase in Histogram['Reflection Lists']:
    30893118        refDict = Histogram['Reflection Lists'][phase]
    … …  
    31153144                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
    31163145        badPeak = False
    3117         for iref,refl in enumerate(refDict['RefList']):
    3118             if 'C' in calcControls[hfx+'histType']:
     3146        # test to see if we are using multiprocessing here
     3147        useMP,ncores = G2mp.InitMP()
     3148        if len(refDict['RefList']) < 100: useMP = False       
     3149        if useMP: # multiprocessing: create a set of initialized Python processes
     3150            MPpool = mp.Pool(ncores,G2mp.InitPwdrProfGlobals,[im,shl,x])
     3151            profArgs = [[] for i in range(ncores)]
     3152        if 'C' in calcControls[hfx+'histType']:
     3153            for iref,refl in enumerate(refDict['RefList']):
    31193154                if im:
    31203155                    h,k,l,m = refl[:4]
    … …  
    31483183                    badPeak = True
    31493184                    continue
    3150                 yc[iBeg:iFin] += refl[11+im]*refl[9+im]*G2pwd.getFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))    #>90% of time spent here
     3185                if useMP:
     3186                    profArgs[iref%ncores].append((refl[5+im],refl,iBeg,iFin,1.))
     3187                else:
     3188                    yc[iBeg:iFin] += refl[11+im]*refl[9+im]*G2pwd.getFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))    #>90% of time spent here
    31513189                if Ka2:
    31523190                    pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
    … …  
    31603198                    elif iBeg > iFin:   #bad peak coeff - skip
    31613199                        continue
    3162                     yc[iBeg:iFin] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))        #and here
    3163             elif 'T' in calcControls[hfx+'histType']:
     3200                    if useMP:
     3201                        profArgs[iref%ncores].append((pos2,refl,iBeg,iFin,kRatio))
     3202                    else:
     3203                        yc[iBeg:iFin] += refl[11+im]*refl[9+im]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))        #and here
     3204        elif 'T' in calcControls[hfx+'histType']:
     3205            for iref,refl in enumerate(refDict['RefList']):
    31643206                h,k,l = refl[:3]
    31653207                Uniq = np.inner(refl[:3],SGMT)
    … …  
    31913233                    badPeak = True
    31923234                    continue
    3193                 yc[iBeg:iFin] += refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))/cw[iBeg:iFin]
     3235                if useMP:
     3236                    profArgs[iref%ncores].append((refl[5+im],refl,iBeg,iFin))
     3237                else:
     3238                    yc[iBeg:iFin] += refl[11+im]*refl[9+im]*G2pwd.getEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))/cw[iBeg:iFin]
    31943239#        print 'profile calc time: %.3fs'%(time.time()-time0)
     3240        if useMP and 'C' in calcControls[hfx+'histType']:
     3241            for y in MPpool.imap_unordered(G2mp.ComputePwdrProfCW,profArgs):
     3242                yc += y
     3243            MPpool.terminate()
     3244        elif useMP:
     3245            for y in MPpool.imap_unordered(G2mp.ComputePwdrProfTOF,profArgs):
     3246                yc += y
     3247            MPpool.terminate()
    31953248    if badPeak:
    31963249        print 'ouch #4 bad profile coefficients yield negative peak width; some reflections skipped'
    3197     #print 'end getPowderProfile t=',time.time()-starttime
     3250    if GSASIIpath.GetConfigValue('debug'):
     3251        print 'getPowderProfile t=',time.time()-starttime
    31983252    return yc,yb
    31993253   
    3200 def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup,dependentVars):
    3201     '''Computes the derivatives of the computed powder pattern with respect to all
    3202     refined parameters
    3203     '''
    3204     #if GSASIIpath.GetConfigValue('debug'):
    3205     #    starttime = time.time()
    3206     #    print 'starting getPowderProfileDerv'
     3254# def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup,dependentVars):
     3255#     '''Computes the derivatives of the computed powder pattern with respect to all
     3256#     refined parameters
     3257#     '''
     3258#     #if GSASIIpath.GetConfigValue('debug'):
     3259#     #    starttime = time.time()
     3260#     #    print 'starting getPowderProfileDerv'
    32073261   
    3208     def cellVaryDerv(pfx,SGData,dpdA):
    3209         if SGData['SGLaue'] in ['-1',]:
    3210             return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
    3211                 [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
    3212         elif SGData['SGLaue'] in ['2/m',]:
    3213             if SGData['SGUniq'] == 'a':
    3214                 return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
    3215             elif SGData['SGUniq'] == 'b':
    3216                 return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
    3217             else:
    3218                 return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
    3219         elif SGData['SGLaue'] in ['mmm',]:
    3220             return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
    3221         elif SGData['SGLaue'] in ['4/m','4/mmm']:
    3222             return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
    3223         elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
    3224             return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
    3225         elif SGData['SGLaue'] in ['3R', '3mR']:
    3226             return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
    3227         elif SGData['SGLaue'] in ['m3m','m3']:
    3228             return [[pfx+'A0',dpdA[0]]]
     3262#     def cellVaryDerv(pfx,SGData,dpdA):
     3263#         if SGData['SGLaue'] in ['-1',]:
     3264#             return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
     3265#                 [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
     3266#         elif SGData['SGLaue'] in ['2/m',]:
     3267#             if SGData['SGUniq'] == 'a':
     3268#                 return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
     3269#             elif SGData['SGUniq'] == 'b':
     3270#                 return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
     3271#             else:
     3272#                 return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
     3273#         elif SGData['SGLaue'] in ['mmm',]:
     3274#             return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
     3275#         elif SGData['SGLaue'] in ['4/m','4/mmm']:
     3276#             return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
     3277#         elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
     3278#             return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
     3279#         elif SGData['SGLaue'] in ['3R', '3mR']:
     3280#             return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
     3281#         elif SGData['SGLaue'] in ['m3m','m3']:
     3282#             return [[pfx+'A0',dpdA[0]]]
    32293283           
    3230     # create a list of dependent variables and set up a dictionary to hold their derivatives
    3231     depDerivDict = {}
    3232     for j in dependentVars:
    3233         depDerivDict[j] = np.zeros(shape=(len(x)))
    3234     #print 'dependent vars',dependentVars
    3235     hId = Histogram['hId']
    3236     hfx = ':%d:'%(hId)
    3237     bakType = calcControls[hfx+'bakType']
    3238     dMdv = np.zeros(shape=(len(varylist),len(x)))
    3239     # do not need dMdv to be a masked array at this point. Moved conversion to later in this routine.
    3240     #dMdv = ma.array(dMdv,mask=np.outer(np.ones(len(varylist)),ma.getmaskarray(x)))      #x is a MaskedArray!
    3241     dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
    3242     if hfx+'Back;0' in varylist: # for now assume that Back;x vars to not appear in constraints
    3243         bBpos = varylist.index(hfx+'Back;0')
    3244         dMdv[bBpos:bBpos+len(dMdb)] += dMdb     #TODO crash if bck parms tossed
    3245     names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
    3246     for name in varylist:
    3247         if 'Debye' in name:
    3248             id = int(name.split(';')[-1])
    3249             parm = name[:int(name.rindex(';'))]
    3250             ip = names.index(parm)
    3251             dMdv[varylist.index(name)] += dMddb[3*id+ip]
    3252     names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
    3253     for name in varylist:
    3254         if 'BkPk' in name:
    3255             parm,id = name.split(';')
    3256             id = int(id)
    3257             if parm in names:
    3258                 ip = names.index(parm)
    3259                 dMdv[varylist.index(name)] += dMdpk[4*id+ip]
    3260     cw = np.diff(ma.getdata(x))
    3261     cw = np.append(cw,cw[-1])
    3262     Ka2 = False #also for TOF!
    3263     if 'C' in calcControls[hfx+'histType']:   
    3264         shl = max(parmDict[hfx+'SH/L'],0.002)
    3265         if hfx+'Lam1' in parmDict.keys():
    3266             wave = parmDict[hfx+'Lam1']
    3267             Ka2 = True
    3268             lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
    3269             kRatio = parmDict[hfx+'I(L2)/I(L1)']
    3270         else:
    3271             wave = parmDict[hfx+'Lam']
    3272     #print '#1 getPowderProfileDerv t=',time.time()-starttime
    3273     for phase in Histogram['Reflection Lists']:
    3274         refDict = Histogram['Reflection Lists'][phase]
    3275         if phase not in Phases:     #skips deleted or renamed phases silently!
    3276             continue
    3277         Phase = Phases[phase]
    3278         SGData = Phase['General']['SGData']
    3279         SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
    3280         im = 0
    3281         if Phase['General'].get('Modulated',False):
    3282             SSGData = Phase['General']['SSGData']
    3283             im = 1  #offset in SS reflection list
    3284             #??
    3285         pId = Phase['pId']
    3286         pfx = '%d::'%(pId)
    3287         phfx = '%d:%d:'%(pId,hId)
    3288         Dij = GetDij(phfx,SGData,parmDict)
    3289         A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
    3290         G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
    3291         GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
    3292         if not Phase['General'].get('doPawley') and not parmDict[phfx+'LeBail']:
    3293             if im:
    3294                 dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
    3295             else:
    3296                 if Phase['General']['Type'] == 'magnetic':
    3297                     dFdvDict = StructureFactorDervMag(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
    3298                 else:
    3299                     dFdvDict = StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
    3300 #            print 'sf-derv time %.3fs'%(time.time()-time0)
    3301             ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
    3302         #print '#2 getPowderProfileDerv t=',time.time()-starttime
    3303         # determine the parameters that will have derivatives computed only at end
    3304         nonatomvarylist = []
    3305         for name in varylist:
    3306             if '::RBV;' not in name:
    3307                 try:
    3308                     aname = name.split(pfx)[1][:2]
    3309                     if aname not in ['Af','dA','AU','RB','AM','Xs','Xc','Ys','Yc','Zs','Zc',    \
    3310                         'Tm','Xm','Ym','Zm','U1','U2','U3']: continue # skip anything not an atom or rigid body param
    3311                 except IndexError:
    3312                     continue
    3313             nonatomvarylist.append(name)
    3314         nonatomdependentVars = []
    3315         for name in dependentVars:
    3316             if '::RBV;' not in name:
    3317                 try:
    3318                     aname = name.split(pfx)[1][:2]
    3319                     if aname not in ['Af','dA','AU','RB','AM','Xs','Xc','Ys','Yc','Zs','Zc',    \
    3320                         'Tm','Xm','Ym','Zm','U1','U2','U3']: continue # skip anything not an atom or rigid body param
    3321                 except IndexError:
    3322                     continue
    3323             nonatomdependentVars.append(name)
    3324         #timelist =  10*[0.0]
    3325         #timestart = 10*[0.0]
    3326         #==========================================================================================
    3327         #==========================================================================================
    3328         for iref,refl in enumerate(refDict['RefList']):
    3329             #timestart[0] = time.time()
    3330             if im:
    3331                 h,k,l,m = refl[:4]
    3332             else:
    3333                 h,k,l = refl[:3]
    3334             Uniq = np.inner(refl[:3],SGMT)
    3335             if 'T' in calcControls[hfx+'histType']:
    3336                 wave = refl[14+im]
    3337             dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,im,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
    3338             if 'C' in calcControls[hfx+'histType']:        #CW powder
    3339                 Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
    3340             else: #'T'OF
    3341                 Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
    3342             iBeg = np.searchsorted(x,refl[5+im]-fmin)
    3343             iFin = np.searchsorted(x,refl[5+im]+fmax)
    3344             if not iBeg+iFin:       #peak below low limit - skip peak
    3345                 continue
    3346             elif not iBeg-iFin:     #peak above high limit - done
    3347                 break
    3348             pos = refl[5+im]
    3349             #itim=0;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
    3350             if 'C' in calcControls[hfx+'histType']:
    3351                 tanth = tand(pos/2.0)
    3352                 costh = cosd(pos/2.0)
    3353                 lenBF = iFin-iBeg
    3354                 dMdpk = np.zeros(shape=(6,lenBF))
    3355                 dMdipk = G2pwd.getdFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))
    3356                 for i in range(5):
    3357                     dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11+im]*refl[9+im]*dMdipk[i]
    3358                 dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
    3359                 if Ka2:
    3360                     pos2 = refl[5+im]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
    3361                     iBeg2 = np.searchsorted(x,pos2-fmin)
    3362                     iFin2 = np.searchsorted(x,pos2+fmax)
    3363                     if iBeg2-iFin2:
    3364                         lenBF2 = iFin2-iBeg2
    3365                         dMdpk2 = np.zeros(shape=(6,lenBF2))
    3366                         dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2]))
    3367                         for i in range(5):
    3368                             dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[11+im]*refl[9+im]*kRatio*dMdipk2[i]
    3369                         dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[11+im]*dMdipk2[0]
    3370                         dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
    3371             else:   #'T'OF
    3372                 lenBF = iFin-iBeg
    3373                 if lenBF < 0:   #bad peak coeff
    3374                     break
    3375                 dMdpk = np.zeros(shape=(6,lenBF))
    3376                 dMdipk = G2pwd.getdEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))
    3377                 for i in range(6):
    3378                     dMdpk[i] += refl[11+im]*refl[9+im]*dMdipk[i]      #cw[iBeg:iFin]*
    3379                 dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}           
    3380             #itim=1;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
    3381             if Phase['General'].get('doPawley'):
    3382                 dMdpw = np.zeros(len(x))
    3383                 try:
    3384                     if im:
    3385                         pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
    3386                     else:
    3387                         pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
    3388                     idx = varylist.index(pIdx)
    3389                     dMdpw[iBeg:iFin] = dervDict['int']/refl[9+im]
    3390                     if Ka2: #not for TOF either
    3391                         dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9+im]
    3392                     dMdv[idx] = dMdpw
    3393                 except: # ValueError:
    3394                     pass
    3395             if 'C' in calcControls[hfx+'histType']:
    3396                 dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY,dpdV = GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict)
    3397                 names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
    3398                     hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
    3399                     hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
    3400                     hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
    3401                     hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
    3402                     hfx+'DisplaceY':[dpdY,'pos'],}
    3403                 if 'Bragg' in calcControls[hfx+'instType']:
    3404                     names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
    3405                         hfx+'SurfRoughB':[dFdAb[1],'int'],})
    3406                 else:
    3407                     names.update({hfx+'Absorption':[dFdAb,'int'],})
    3408             else:   #'T'OF
    3409                 dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV = GetReflPosDerv(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)
    3410                 names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
    3411                     hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
    3412                     hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4+im],'gam'],hfx+'Y':[refl[4+im]**2,'gam'],
    3413                     hfx+'alpha':[1./refl[4+im],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4+im]**4,'bet'],
    3414                     hfx+'beta-q':[1./refl[4+im]**2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4+im]**2,'sig'],
    3415                     hfx+'sig-2':[refl[4+im]**4,'sig'],hfx+'sig-q':[1./refl[4+im]**2,'sig'],
    3416                     hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
    3417             #itim=2;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
    3418             for name in names:
    3419                 item = names[name]
    3420                 if name in varylist:
    3421                     dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
    3422                     if Ka2 and iFin2-iBeg2:
    3423                         dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
    3424                 elif name in dependentVars:
    3425                     depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
    3426                     if Ka2 and iFin2-iBeg2:
    3427                         depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
    3428             for iPO in dIdPO:
    3429                 if iPO in varylist:
    3430                     dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
    3431                     if Ka2 and iFin2-iBeg2:
    3432                         dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
    3433                 elif iPO in dependentVars:
    3434                     depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
    3435                     if Ka2 and iFin2-iBeg2:
    3436                         depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
    3437             #itim=3;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
    3438             for i,name in enumerate(['omega','chi','phi']):
    3439                 aname = pfx+'SH '+name
    3440                 if aname in varylist:
    3441                     dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
    3442                     if Ka2 and iFin2-iBeg2:
    3443                         dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
    3444                 elif aname in dependentVars:
    3445                     depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
    3446                     if Ka2 and iFin2-iBeg2:
    3447                         depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
    3448             for iSH in dFdODF:
    3449                 if iSH in varylist:
    3450                     dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
    3451                     if Ka2 and iFin2-iBeg2:
    3452                         dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
    3453                 elif iSH in dependentVars:
    3454                     depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
    3455                     if Ka2 and iFin2-iBeg2:
    3456                         depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
    3457             cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
    3458             for name,dpdA in cellDervNames:
    3459                 if name in varylist:
    3460                     dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
    3461                     if Ka2 and iFin2-iBeg2:
    3462                         dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
    3463                 elif name in dependentVars: #need to scale for mixed phase constraints?
    3464                     depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
    3465                     if Ka2 and iFin2-iBeg2:
    3466                         depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
    3467             dDijDict = GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict)
    3468             for name in dDijDict:
    3469                 if name in varylist:
    3470                     dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
    3471                     if Ka2 and iFin2-iBeg2:
    3472                         dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
    3473                 elif name in dependentVars:
    3474                     depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
    3475                     if Ka2 and iFin2-iBeg2:
    3476                         depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
    3477             #itim=4;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
    3478             for i,name in enumerate([pfx+'mV0',pfx+'mV1',pfx+'mV2']):
    3479                 if name in varylist:
    3480                     dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict['pos']
    3481                     if Ka2 and iFin2-iBeg2:
    3482                         dMdv[varylist.index(name)][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
    3483                 elif name in dependentVars:
    3484                     depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict['pos']
    3485                     if Ka2 and iFin2-iBeg2:
    3486                         depDerivDict[name][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
    3487             if 'C' in calcControls[hfx+'histType']:
    3488                 sigDict,gamDict = GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
    3489             else:   #'T'OF
    3490                 sigDict,gamDict = GetSampleSigGamDerv(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
    3491             for name in gamDict:
    3492                 if name in varylist:
    3493                     dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
    3494                     if Ka2 and iFin2-iBeg2:
    3495                         dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
    3496                 elif name in dependentVars:
    3497                     depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
    3498                     if Ka2 and iFin2-iBeg2:
    3499                         depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
    3500             for name in sigDict:
    3501                 if name in varylist:
    3502                     dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
    3503                     if Ka2 and iFin2-iBeg2:
    3504                         dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
    3505                 elif name in dependentVars:
    3506                     depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
    3507                     if Ka2 and iFin2-iBeg2:
    3508                         depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
    3509             for name in ['BabA','BabU']:
    3510                 if refl[9+im]:
    3511                     if phfx+name in varylist:
    3512                         dMdv[varylist.index(phfx+name)][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict['int']/refl[9+im]
    3513                         if Ka2 and iFin2-iBeg2:
    3514                             dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict2['int']/refl[9+im]
    3515                     elif phfx+name in dependentVars:                   
    3516                         depDerivDict[phfx+name][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict['int']/refl[9+im]
    3517                         if Ka2 and iFin2-iBeg2:
    3518                             depDerivDict[phfx+name][iBeg2:iFin2] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict2['int']/refl[9+im]                 
    3519             #itim=5;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
    3520             if not Phase['General'].get('doPawley') and not parmDict[phfx+'LeBail']:
    3521                 #do atom derivatives -  for RB,F,X & U so far - how do I scale mixed phase constraints?
    3522                 corr = 0.
    3523                 corr2 = 0.
    3524                 if refl[9+im]:             
    3525                     corr = dervDict['int']/refl[9+im]
    3526                     #if Ka2 and iFin2-iBeg2:
    3527                     #    corr2 = dervDict2['int']/refl[9+im]
    3528                 #itim=6;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
    3529                 for name in nonatomvarylist:
    3530                     dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
    3531                     if Ka2 and iFin2-iBeg2:
    3532                        dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
    3533                 #itim=7;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
    3534                 for name in nonatomdependentVars:
    3535                    depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
    3536                    if Ka2 and iFin2-iBeg2:
    3537                        depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
    3538                 #itim=8;timelist[itim] += time.time()-timestart[itim]   
    3539     #        print 'profile derv time: %.3fs'%(time.time()-time0)
    3540     # now process derivatives in constraints
    3541     #print '#3 getPowderProfileDerv t=',time.time()-starttime
    3542     #print timelist,sum(timelist)
    3543     dMdv[:,ma.getmaskarray(x)] = 0.  # instead of masking, zero out masked values
    3544     dMdv = ma.array(dMdv,mask=np.outer(np.ones(len(varylist)),ma.getmaskarray(x)))      #x is a MaskedArray!
    3545     G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
    3546     #if GSASIIpath.GetConfigValue('debug'):
    3547     #    print 'end getPowderProfileDerv t=',time.time()-starttime
    3548     return dMdv
     3284#     # create a list of dependent variables and set up a dictionary to hold their derivatives
     3285#     depDerivDict = {}
     3286#     for j in dependentVars:
     3287#         depDerivDict[j] = np.zeros(shape=(len(x)))
     3288#     #print 'dependent vars',dependentVars
     3289#     hId = Histogram['hId']
     3290#     hfx = ':%d:'%(hId)
     3291#     bakType = calcControls[hfx+'bakType']
     3292#     dMdv = np.zeros(shape=(len(varylist),len(x)))
     3293#     # do not need dMdv to be a masked array at this point. Moved conversion to later in this routine.
     3294#     #dMdv = ma.array(dMdv,mask=np.outer(np.ones(len(varylist)),ma.getmaskarray(x)))      #x is a MaskedArray!
     3295#     dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,calcControls[hfx+'histType'],x)
     3296#     if hfx+'Back;0' in varylist: # for now assume that Back;x vars to not appear in constraints
     3297#         bBpos = varylist.index(hfx+'Back;0')
     3298#         dMdv[bBpos:bBpos+len(dMdb)] += dMdb     #TODO crash if bck parms tossed
     3299#     names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
     3300#     for name in varylist:
     3301#         if 'Debye' in name:
     3302#             id = int(name.split(';')[-1])
     3303#             parm = name[:int(name.rindex(';'))]
     3304#             ip = names.index(parm)
     3305#             dMdv[varylist.index(name)] += dMddb[3*id+ip]
     3306#     names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
     3307#     for name in varylist:
     3308#         if 'BkPk' in name:
     3309#             parm,id = name.split(';')
     3310#             id = int(id)
     3311#             if parm in names:
     3312#                 ip = names.index(parm)
     3313#                 dMdv[varylist.index(name)] += dMdpk[4*id+ip]
     3314#     cw = np.diff(ma.getdata(x))
     3315#     cw = np.append(cw,cw[-1])
     3316#     Ka2 = False #also for TOF!
     3317#     if 'C' in calcControls[hfx+'histType']:   
     3318#         shl = max(parmDict[hfx+'SH/L'],0.002)
     3319#         if hfx+'Lam1' in parmDict.keys():
     3320#             wave = parmDict[hfx+'Lam1']
     3321#             Ka2 = True
     3322#             lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
     3323#             kRatio = parmDict[hfx+'I(L2)/I(L1)']
     3324#         else:
     3325#             wave = parmDict[hfx+'Lam']
     3326#     #print '#1 getPowderProfileDerv t=',time.time()-starttime
     3327#     for phase in Histogram['Reflection Lists']:
     3328#         refDict = Histogram['Reflection Lists'][phase]
     3329#         if phase not in Phases:     #skips deleted or renamed phases silently!
     3330#             continue
     3331#         Phase = Phases[phase]
     3332#         SGData = Phase['General']['SGData']
     3333#         SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     3334#         im = 0
     3335#         if Phase['General'].get('Modulated',False):
     3336#             SSGData = Phase['General']['SSGData']
     3337#             im = 1  #offset in SS reflection list
     3338#             #??
     3339#         pId = Phase['pId']
     3340#         pfx = '%d::'%(pId)
     3341#         phfx = '%d:%d:'%(pId,hId)
     3342#         Dij = GetDij(phfx,SGData,parmDict)
     3343#         A = [parmDict[pfx+'A%d'%(i)]+Dij[i] for i in range(6)]
     3344#         G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
     3345#         GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
     3346#         if not Phase['General'].get('doPawley') and not parmDict[phfx+'LeBail']:
     3347#             if im:
     3348#                 dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
     3349#             else:
     3350#                 if Phase['General']['Type'] == 'magnetic':
     3351#                     dFdvDict = StructureFactorDervMag(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
     3352#                 else:
     3353#                     dFdvDict = StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
     3354# #            print 'sf-derv time %.3fs'%(time.time()-time0)
     3355#             ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
     3356#         #print '#2 getPowderProfileDerv t=',time.time()-starttime
     3357#         # determine the parameters that will have derivatives computed only at end
     3358#         nonatomvarylist = []
     3359#         for name in varylist:
     3360#             if '::RBV;' not in name:
     3361#                 try:
     3362#                     aname = name.split(pfx)[1][:2]
     3363#                     if aname not in ['Af','dA','AU','RB','AM','Xs','Xc','Ys','Yc','Zs','Zc',    \
     3364#                         'Tm','Xm','Ym','Zm','U1','U2','U3']: continue # skip anything not an atom or rigid body param
     3365#                 except IndexError:
     3366#                     continue
     3367#             nonatomvarylist.append(name)
     3368#         nonatomdependentVars = []
     3369#         for name in dependentVars:
     3370#             if '::RBV;' not in name:
     3371#                 try:
     3372#                     aname = name.split(pfx)[1][:2]
     3373#                     if aname not in ['Af','dA','AU','RB','AM','Xs','Xc','Ys','Yc','Zs','Zc',    \
     3374#                         'Tm','Xm','Ym','Zm','U1','U2','U3']: continue # skip anything not an atom or rigid body param
     3375#                 except IndexError:
     3376#                     continue
     3377#             nonatomdependentVars.append(name)
     3378#         #==========================================================================================
     3379#         #==========================================================================================
     3380#         for iref,refl in enumerate(refDict['RefList']):
     3381#             if im:
     3382#                 h,k,l,m = refl[:4]
     3383#             else:
     3384#                 h,k,l = refl[:3]
     3385#             Uniq = np.inner(refl[:3],SGMT)
     3386#             if 'T' in calcControls[hfx+'histType']:
     3387#                 wave = refl[14+im]
     3388#             dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb,dFdEx = GetIntensityDerv(refl,im,wave,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
     3389#             if 'C' in calcControls[hfx+'histType']:        #CW powder
     3390#                 Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5+im],refl[6+im],refl[7+im],shl)
     3391#             else: #'T'OF
     3392#                 Wd,fmin,fmax = G2pwd.getWidthsTOF(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im])
     3393#             iBeg = np.searchsorted(x,refl[5+im]-fmin)
     3394#             iFin = np.searchsorted(x,refl[5+im]+fmax)
     3395#             if not iBeg+iFin:       #peak below low limit - skip peak
     3396#                 continue
     3397#             elif not iBeg-iFin:     #peak above high limit - done
     3398#                 break
     3399#             pos = refl[5+im]
     3400#             #itim=0;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
     3401#             if 'C' in calcControls[hfx+'histType']:
     3402#                 tanth = tand(pos/2.0)
     3403#                 costh = cosd(pos/2.0)
     3404#                 lenBF = iFin-iBeg
     3405#                 dMdpk = np.zeros(shape=(6,lenBF))
     3406#                 dMdipk = G2pwd.getdFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))
     3407#                 for i in range(5):
     3408#                     dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11+im]*refl[9+im]*dMdipk[i]
     3409#                 dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
     3410#                 if Ka2:
     3411#                     pos2 = refl[5+im]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
     3412#                     iBeg2 = np.searchsorted(x,pos2-fmin)
     3413#                     iFin2 = np.searchsorted(x,pos2+fmax)
     3414#                     if iBeg2-iFin2:
     3415#                         lenBF2 = iFin2-iBeg2
     3416#                         dMdpk2 = np.zeros(shape=(6,lenBF2))
     3417#                         dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg2:iFin2]))
     3418#                         for i in range(5):
     3419#                             dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[11+im]*refl[9+im]*kRatio*dMdipk2[i]
     3420#                         dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[11+im]*dMdipk2[0]
     3421#                         dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
     3422#             else:   #'T'OF
     3423#                 lenBF = iFin-iBeg
     3424#                 if lenBF < 0:   #bad peak coeff
     3425#                     break
     3426#                 dMdpk = np.zeros(shape=(6,lenBF))
     3427#                 dMdipk = G2pwd.getdEpsVoigt(refl[5+im],refl[12+im],refl[13+im],refl[6+im],refl[7+im],ma.getdata(x[iBeg:iFin]))
     3428#                 for i in range(6):
     3429#                     dMdpk[i] += refl[11+im]*refl[9+im]*dMdipk[i]      #cw[iBeg:iFin]*
     3430#                 dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]}           
     3431#             #itim=1;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
     3432#             if Phase['General'].get('doPawley'):
     3433#                 dMdpw = np.zeros(len(x))
     3434#                 try:
     3435#                     if im:
     3436#                         pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d,%d'%(h,k,l,m)])
     3437#                     else:
     3438#                         pIdx = pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
     3439#                     idx = varylist.index(pIdx)
     3440#                     dMdpw[iBeg:iFin] = dervDict['int']/refl[9+im]
     3441#                     if Ka2: #not for TOF either
     3442#                         dMdpw[iBeg2:iFin2] += dervDict2['int']/refl[9+im]
     3443#                     dMdv[idx] = dMdpw
     3444#                 except: # ValueError:
     3445#                     pass
     3446#             if 'C' in calcControls[hfx+'histType']:
     3447#                 dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY,dpdV = GetReflPosDerv(refl,im,wave,A,pfx,hfx,calcControls,parmDict)
     3448#                 names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
     3449#                     hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
     3450#                     hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
     3451#                     hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
     3452#                     hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
     3453#                     hfx+'DisplaceY':[dpdY,'pos'],}
     3454#                 if 'Bragg' in calcControls[hfx+'instType']:
     3455#                     names.update({hfx+'SurfRoughA':[dFdAb[0],'int'],
     3456#                         hfx+'SurfRoughB':[dFdAb[1],'int'],})
     3457#                 else:
     3458#                     names.update({hfx+'Absorption':[dFdAb,'int'],})
     3459#             else:   #'T'OF
     3460#                 dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV = GetReflPosDerv(refl,im,0.0,A,pfx,hfx,calcControls,parmDict)
     3461#                 names = {hfx+'Scale':[dIdsh,'int'],phfx+'Scale':[dIdsp,'int'],
     3462#                     hfx+'difC':[dpdDC,'pos'],hfx+'difA':[dpdDA,'pos'],hfx+'difB':[dpdDB,'pos'],
     3463#                     hfx+'Zero':[dpdZ,'pos'],hfx+'X':[refl[4+im],'gam'],hfx+'Y':[refl[4+im]**2,'gam'],
     3464#                     hfx+'alpha':[1./refl[4+im],'alp'],hfx+'beta-0':[1.0,'bet'],hfx+'beta-1':[1./refl[4+im]**4,'bet'],
     3465#                     hfx+'beta-q':[1./refl[4+im]**2,'bet'],hfx+'sig-0':[1.0,'sig'],hfx+'sig-1':[refl[4+im]**2,'sig'],
     3466#                     hfx+'sig-2':[refl[4+im]**4,'sig'],hfx+'sig-q':[1./refl[4+im]**2,'sig'],
     3467#                     hfx+'Absorption':[dFdAb,'int'],phfx+'Extinction':[dFdEx,'int'],}
     3468#             #itim=2;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
     3469#             for name in names:
     3470#                 item = names[name]
     3471#                 if name in varylist:
     3472#                     dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
     3473#                     if Ka2 and iFin2-iBeg2:
     3474#                         dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
     3475#                 elif name in dependentVars:
     3476#                     depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
     3477#                     if Ka2 and iFin2-iBeg2:
     3478#                         depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
     3479#             for iPO in dIdPO:
     3480#                 if iPO in varylist:
     3481#                     dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
     3482#                     if Ka2 and iFin2-iBeg2:
     3483#                         dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
     3484#                 elif iPO in dependentVars:
     3485#                     depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
     3486#                     if Ka2 and iFin2-iBeg2:
     3487#                         depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
     3488#             #itim=3;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
     3489#             for i,name in enumerate(['omega','chi','phi']):
     3490#                 aname = pfx+'SH '+name
     3491#                 if aname in varylist:
     3492#                     dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
     3493#                     if Ka2 and iFin2-iBeg2:
     3494#                         dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
     3495#                 elif aname in dependentVars:
     3496#                     depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
     3497#                     if Ka2 and iFin2-iBeg2:
     3498#                         depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
     3499#             for iSH in dFdODF:
     3500#                 if iSH in varylist:
     3501#                     dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
     3502#                     if Ka2 and iFin2-iBeg2:
     3503#                         dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
     3504#                 elif iSH in dependentVars:
     3505#                     depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
     3506#                     if Ka2 and iFin2-iBeg2:
     3507#                         depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
     3508#             cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
     3509#             for name,dpdA in cellDervNames:
     3510#                 if name in varylist:
     3511#                     dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
     3512#                     if Ka2 and iFin2-iBeg2:
     3513#                         dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
     3514#                 elif name in dependentVars: #need to scale for mixed phase constraints?
     3515#                     depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
     3516#                     if Ka2 and iFin2-iBeg2:
     3517#                         depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
     3518#             dDijDict = GetHStrainShiftDerv(refl,im,SGData,phfx,hfx,calcControls,parmDict)
     3519#             for name in dDijDict:
     3520#                 if name in varylist:
     3521#                     dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
     3522#                     if Ka2 and iFin2-iBeg2:
     3523#                         dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
     3524#                 elif name in dependentVars:
     3525#                     depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
     3526#                     if Ka2 and iFin2-iBeg2:
     3527#                         depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
     3528#             #itim=4;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
     3529#             for i,name in enumerate([pfx+'mV0',pfx+'mV1',pfx+'mV2']):
     3530#                 if name in varylist:
     3531#                     dMdv[varylist.index(name)][iBeg:iFin] += dpdV[i]*dervDict['pos']
     3532#                     if Ka2 and iFin2-iBeg2:
     3533#                         dMdv[varylist.index(name)][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
     3534#                 elif name in dependentVars:
     3535#                     depDerivDict[name][iBeg:iFin] += dpdV[i]*dervDict['pos']
     3536#                     if Ka2 and iFin2-iBeg2:
     3537#                         depDerivDict[name][iBeg2:iFin2] += dpdV[i]*dervDict2['pos']
     3538#             if 'C' in calcControls[hfx+'histType']:
     3539#                 sigDict,gamDict = GetSampleSigGamDerv(refl,im,wave,G,GB,SGData,hfx,phfx,calcControls,parmDict)
     3540#             else:   #'T'OF
     3541#                 sigDict,gamDict = GetSampleSigGamDerv(refl,im,0.0,G,GB,SGData,hfx,phfx,calcControls,parmDict)
     3542#             for name in gamDict:
     3543#                 if name in varylist:
     3544#                     dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
     3545#                     if Ka2 and iFin2-iBeg2:
     3546#                         dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
     3547#                 elif name in dependentVars:
     3548#                     depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
     3549#                     if Ka2 and iFin2-iBeg2:
     3550#                         depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
     3551#             for name in sigDict:
     3552#                 if name in varylist:
     3553#                     dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
     3554#                     if Ka2 and iFin2-iBeg2:
     3555#                         dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
     3556#                 elif name in dependentVars:
     3557#                     depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
     3558#                     if Ka2 and iFin2-iBeg2:
     3559#                         depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
     3560#             for name in ['BabA','BabU']:
     3561#                 if refl[9+im]:
     3562#                     if phfx+name in varylist:
     3563#                         dMdv[varylist.index(phfx+name)][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict['int']/refl[9+im]
     3564#                         if Ka2 and iFin2-iBeg2:
     3565#                             dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict2['int']/refl[9+im]
     3566#                     elif phfx+name in dependentVars:                   
     3567#                         depDerivDict[phfx+name][iBeg:iFin] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict['int']/refl[9+im]
     3568#                         if Ka2 and iFin2-iBeg2:
     3569#                             depDerivDict[phfx+name][iBeg2:iFin2] += parmDict[phfx+'Scale']*dFdvDict[phfx+name][iref]*dervDict2['int']/refl[9+im]                 
     3570#             #itim=5;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
     3571#             if not Phase['General'].get('doPawley') and not parmDict[phfx+'LeBail']:
     3572#                 #do atom derivatives -  for RB,F,X & U so far - how do I scale mixed phase constraints?
     3573#                 corr = 0.
     3574#                 corr2 = 0.
     3575#                 if refl[9+im]:             
     3576#                     corr = dervDict['int']/refl[9+im]
     3577#                     #if Ka2 and iFin2-iBeg2:
     3578#                     #    corr2 = dervDict2['int']/refl[9+im]
     3579#                 #itim=6;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
     3580#                 for name in nonatomvarylist:
     3581#                     dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
     3582#                     if Ka2 and iFin2-iBeg2:
     3583#                        dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
     3584#                 #itim=7;timelist[itim] += time.time()-timestart[itim]; timestart[itim+1] = time.time()
     3585#                 for name in nonatomdependentVars:
     3586#                    depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
     3587#                    if Ka2 and iFin2-iBeg2:
     3588#                        depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
     3589#                 #itim=8;timelist[itim] += time.time()-timestart[itim]   
     3590#     #        print 'profile derv time: %.3fs'%(time.time()-time0)
     3591#     # now process derivatives in constraints
     3592#     #print '#3 getPowderProfileDerv t=',time.time()-starttime
     3593#     #print timelist,sum(timelist)
     3594#     dMdv[:,ma.getmaskarray(x)] = 0.  # instead of masking, zero out masked values
     3595#     dMdv = ma.array(dMdv,mask=np.outer(np.ones(len(varylist)),ma.getmaskarray(x)))      #x is a MaskedArray!
     3596#     G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
     3597#     #if GSASIIpath.GetConfigValue('debug'):
     3598#     #    print 'end getPowderProfileDerv t=',time.time()-starttime
     3599#     return dMdv
    35493600
    35503601def getPowderProfileDervMP(args):
    … …  
    40784129            xB = np.searchsorted(x,Limits[0])
    40794130            xF = np.searchsorted(x,Limits[1])+1
    4080             ######################################################################
    4081             #import GSASIImpsubs as G2mp
    4082             #G2mp.InitMP()
    4083             useMP = False
    4084             ncores = GSASIIpath.GetConfigValue('Multiprocessing_cores')
    4085             #useMP = G2mp.useMP #and len(refDict['RefList']) > 100
    4086             if GSASIIpath.GetConfigValue('debug'):
    4087                 starttime = time.time()
    4088 #                print 'starting getPowderProfileDerv'
    4089             #useMP = True
    4090             if useMP and ncores > 1:
    4091                 import multiprocessing as mp
    4092                 print 'mp with ',ncores,'cores'
     4131            useMP,ncores = G2mp.InitMP()
     4132            if GSASIIpath.GetConfigValue('debug'): starttime = time.time()
     4133            if useMP:
    40934134                MPpool = mp.Pool(ncores)
    40944135                dMdvh = None
    … …  
    41024143                       dMdvh += dmdv
    41034144            else:
    4104                 #dMdvh = getPowderProfileDervMP([parmDict,x[xB:xF],
    4105                 #    varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup])
    4106                 dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
    4107                     varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup,dependentVars)
    4108             if GSASIIpath.GetConfigValue('debug'):
    4109                 print 'end getPowderProfileDerv t=',time.time()-starttime
    4110             #import cPickle
    4111             #fp = open('/tmp/hess.pkl','w')
    4112             #cPickle.dump(dMdvh,fp,1)
    4113             #fp.close()
    4114             ######################################################################
     4145                dMdvh = getPowderProfileDervMP([parmDict,x[xB:xF],
     4146                    varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup,dependentVars])
     4147                #dMdvh = getPowderProfileDerv(parmDict,x[xB:xF],
     4148                #    varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup,dependentVars)
     4149            if GSASIIpath.GetConfigValue('debug'): print 'getPowderProfileDerv t=',time.time()-starttime
    41154150            Wt = ma.sqrt(W[xB:xF])[nxs,:]
    41164151            Dy = dy[xB:xF][nxs,:]
    … …  
    43874422        M = np.concatenate((M,np.sqrt(pWt)*pVals))
    43884423    return M
    4389                        

  • trunk/config_example.py

    r3033 r3041  
    153153'''
    154154
    155 Multiprocessing_cores = -1
     155Multiprocessing_cores = 0
    156156''' Specifies the number of cores to use when performing multicore computing. A number less
    157 than zero (default) causes the recommended number of cores [using multiprocessing.cpu_count()]
     157than zero causes the recommended number of cores [using multiprocessing.cpu_count()/2]
    158158to be used. Setting this number to 0 or 1 avoids use of the multiprocessing module: all
    159159computations are performed in-line.
Note: See TracChangeset for help on using the changeset viewer.