Changeset 512


Ignore:
Timestamp:
Mar 5, 2012 12:34:16 PM (10 years ago)
Author:
toby
Message:

implement reflection-grouped multiprocessing

File:
1 edited

Legend:

Unmodified
Added
Removed
  • MPbranch/GSASIIstruct.py

    r507 r512  
    4141    '''
    4242    Controls = {'deriv type':'analytic Hessian','max cyc':3,
    43                 'max Hprocess':2,
    44                 'max Rprocess':1,
     43                'max Hprocess':1,
     44                'max Rprocess':2,
    4545                'min dM/M':0.0001,'shift factor':1.}
    4646    file = open(GPXfile,'rb')
     
    20002000    if pfx+'SHorder' in parmDict:
    20012001        Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
    2002     refl[13] = Icorr       
     2002#    refl[13] = Icorr       
     2003    return Icorr       
    20032004   
    20042005def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
     
    22762277                Histogram[phfx+'Nref'] = len(refList)
    22772278               
    2278 def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
    2279    
     2279def ComputeReflectionProfile(arg):
    22802280    def GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
    22812281        U = parmDict[hfx+'U']
     
    22902290        gam = max(0.001,gam)
    22912291        return sig,gam
    2292                
     2292
     2293    (RefNum,refList,x,Phase,calcControls,wave,G,g,GB,Vst,
     2294     parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio) = arg
     2295    yc = np.zeros_like(x)
     2296    iBegO = None
     2297    SGData = Phase['General']['SGData']
     2298    shl = max(parmDict[hfx+'SH/L'],0.002)
     2299    if 'C' in calcControls[hfx+'histType']:
     2300        for refl in refList:
     2301            refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
     2302            Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
     2303            refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
     2304            refl[6:8] = GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict)    #peak sig & gam
     2305            refl[13] = Vst*Lorenz * GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)    #puts corrections in refl[13]
     2306            Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
     2307            iBeg = np.searchsorted(x,refl[5]-fmin)
     2308            iFin = np.searchsorted(x,refl[5]+fmax)
     2309            if not iBeg+iFin:       #peak below low limit - skip peak
     2310                continue
     2311            elif not iBeg-iFin:     #peak above high limit - done
     2312                break
     2313            yc[iBeg:iFin] += refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
     2314            if iBegO is None:
     2315                iBegO = iBeg
     2316                iFinO = iFin
     2317            else:
     2318                iBegO = min(iBegO,iBeg)
     2319                iFinO = max(iFinO,iFin)
     2320            if Ka2:
     2321                pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
     2322                Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
     2323                iBeg2 = np.searchsorted(x,pos2-fmin)
     2324                iFin2 = np.searchsorted(x,pos2+fmax)
     2325                if not iBeg2+iFin2:       #peak below low limit - skip peak
     2326                    continue
     2327                elif not iBeg2-iFin2:     #peak above high limit - done
     2328                    continue
     2329                yc[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
     2330                iBegO = min(iBegO,iBeg2)
     2331                iFinO = max(iFinO,iFin2)
     2332    elif 'T' in calcControls[hfx+'histType']:
     2333        print 'TOF Undefined at present'
     2334        raise Exception    #no TOF yet
     2335
     2336    if iBegO is None: return None # no reflections were computed
     2337    return RefNum,refList,iBegO,iFinO,yc[iBegO:iFinO]
     2338
     2339def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):               
     2340    mpPool = mp.Pool()
    22932341    hId = Histogram['hId']
    22942342    hfx = ':%d:'%(hId)
     
    22992347    if 'C' in calcControls[hfx+'histType']:   
    23002348        shl = max(parmDict[hfx+'SH/L'],0.002)
    2301         Ka2 = False
    23022349        if hfx+'Lam1' in parmDict.keys():
    23032350            wave = parmDict[hfx+'Lam1']
     
    23072354        else:
    23082355            wave = parmDict[hfx+'Lam']
     2356            Ka2 = False
     2357            lamRatio = None
     2358            kRatio = None
    23092359    else:
    23102360        print 'TOF Undefined at present'
     
    23242374        if 'Pawley' not in Phase['General']['Type']:
    23252375            refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
    2326         for refl in refList:
    2327             if 'C' in calcControls[hfx+'histType']:
     2376        else:
     2377            for refl in refList:
    23282378                h,k,l = refl[:3]
    2329                 refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
    2330                 Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
    2331                 refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
    2332                 refl[6:8] = GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict)    #peak sig & gam
    2333                 GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)    #puts corrections in refl[13]
    2334                 refl[13] *= Vst*Lorenz
    2335                 if 'Pawley' in Phase['General']['Type']:
    2336                     try:
    2337                         refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
    2338                     except KeyError:
    2339 #                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
    2340                         continue
    2341                 Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
    2342                 iBeg = np.searchsorted(x,refl[5]-fmin)
    2343                 iFin = np.searchsorted(x,refl[5]+fmax)
    2344                 if not iBeg+iFin:       #peak below low limit - skip peak
     2379                try:
     2380                    refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
     2381                except KeyError:
     2382                    #print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
    23452383                    continue
    2346                 elif not iBeg-iFin:     #peak above high limit - done
    2347                     break
    2348                 yc[iBeg:iFin] += refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
    2349                 if Ka2:
    2350                     pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
    2351                     Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
    2352                     iBeg = np.searchsorted(x,pos2-fmin)
    2353                     iFin = np.searchsorted(x,pos2+fmax)
    2354                     if not iBeg+iFin:       #peak below low limit - skip peak
    2355                         continue
    2356                     elif not iBeg-iFin:     #peak above high limit - done
    2357                         return yc,yb
    2358                     yc[iBeg:iFin] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
    2359             elif 'T' in calcControls[hfx+'histType']:
    2360                 print 'TOF Undefined at present'
    2361                 raise Exception    #no TOF yet
     2384
     2385        chunk = 200
     2386        if (calcControls['max Hprocess'] > 1 or
     2387            calcControls['max Rprocess'] == 1 or
     2388            mp.cpu_count() < 2 or len(refList) <  2 * chunk):
     2389            print 'single thread function',hId,len(refList)
     2390            res = ComputeReflectionProfile(
     2391                (0,refList,x,Phase,calcControls,wave,G,g,GB,Vst,
     2392                 parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio))
     2393            if not res: continue
     2394            RefNum,refls,iBeg,iFin,yct = res
     2395            yc[iBeg:iFin] += yct
     2396        else:
     2397            print 'multiprocess function',hId,len(refList)
     2398            argList = []
     2399            for i in range(0,len(refList),chunk):
     2400                argList.append(
     2401                    (i,refList[i:i+chunk],x,Phase,calcControls,wave,G,g,GB,Vst,
     2402                     parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio))
     2403            for res in mpPool.imap_unordered(ComputeReflectionProfile,argList):
     2404                if not res: continue
     2405                RefNum,refls,iBeg,iFin,yct = res
     2406                yc[iBeg:iFin] += yct
     2407                for i,refl in enumerate(refls):
     2408                    refList[RefNum+i] = refl
     2409            mpPool.close()
     2410            mpPool.join()
    23622411    return yc,yb
    2363    
    2364 def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
    2365    
     2412
     2413def ComputeReflectionDerivative(arg):
    23662414    def cellVaryDerv(pfx,SGData,dpdA):
    23672415        if SGData['SGLaue'] in ['-1',]:
     
    23882436#            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]]
    23892437            return [[pfx+'A0',dpdA[0]]]
     2438
     2439    (RefNum,refList,x,dependentVars,varylist,dFdvDict,
     2440     Phase,calcControls,wave,G,g,GB,A,
     2441     parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio) = arg
     2442    dMdv = np.zeros(shape=(len(varylist),len(x)))
     2443    depDerivDict = {}
     2444    for j in dependentVars:
     2445        depDerivDict[j] = np.zeros(shape=(len(x)))
     2446    SGData = Phase['General']['SGData']
     2447    shl = max(parmDict[hfx+'SH/L'],0.002)
     2448    dx = x[1]-x[0]
     2449    iBegO = None
     2450    for irefS,refl in enumerate(refList):
     2451        iref = irefS + RefNum       
     2452        if 'C' in calcControls[hfx+'histType']:        #CW powder
     2453            dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
     2454            Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
     2455            iBeg = np.searchsorted(x,refl[5]-fmin)
     2456            iFin = np.searchsorted(x,refl[5]+fmax)
     2457            if iBegO is None:
     2458                iBegO = iBeg
     2459                iFinO = iFin
     2460            else:
     2461                iBegO = min(iBegO,iBeg)
     2462                iFinO = max(iFinO,iFin)
     2463            if not iBeg+iFin:       #peak below low limit - skip peak
     2464                continue
     2465                #return None
     2466            elif not iBeg-iFin:     #peak above high limit - done
     2467                break
     2468                #return None
     2469            pos = refl[5]
     2470            tanth = tand(pos/2.0)
     2471            costh = cosd(pos/2.0)
     2472            lenBF = iFin-iBeg
     2473            dMdpk = np.zeros(shape=(6,lenBF))
     2474            dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
     2475            for i in range(1,5):
     2476                dMdpk[i] += 100.*dx*refl[13]*refl[9]*dMdipk[i]
     2477            dMdpk[0] += 100.*dx*refl[13]*refl[9]*dMdipk[0]
     2478            dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
     2479            if Ka2:
     2480                pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
     2481                kdelt = int((pos2-refl[5])/dx)               
     2482                iBeg2 = min(lenX,iBeg+kdelt)
     2483                iFin2 = min(lenX,iFin+kdelt)
     2484                if iBeg2-iFin2:
     2485                    iBegO = min(iBegO,iBeg2)
     2486                    iFinO = max(iFinO,iFin2)
     2487                    lenBF2 = iFin2-iBeg2
     2488                    dMdpk2 = np.zeros(shape=(6,lenBF2))
     2489                    dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
     2490                    for i in range(1,5):
     2491                        dMdpk2[i] = 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[i]
     2492                    dMdpk2[0] = 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[0]
     2493                    dMdpk2[5] = 100.*dx*refl[13]*dMdipk2[0]
     2494                    dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
     2495            if 'Pawley' in Phase['General']['Type']:
     2496                    try:
     2497                        idx = varylist.index(pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%refl[:3]]))
     2498                        dMdv[idx][iBeg:iFin] = dervDict['int']/refl[9]
     2499                        if Ka2:
     2500                            dMdv[idx][iBeg2:iFin2] = dervDict2['int']/refl[9]
     2501                        # Assuming Pawley variables not in constraints
     2502                    except ValueError:
     2503                        pass
     2504            dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
     2505            names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
     2506                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
     2507                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
     2508                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
     2509                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
     2510                    hfx+'DisplaceY':[dpdY,'pos'],}
     2511            for name in names:
     2512                item = names[name]
     2513                if name in varylist:
     2514                    dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
     2515                    if Ka2:
     2516                        dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
     2517                elif name in dependentVars:
     2518                    if Ka2:
     2519                        depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
     2520                    depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
     2521            for iPO in dIdPO:
     2522                if iPO in varylist:
     2523                    dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
     2524                    if Ka2:
     2525                        dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
     2526                elif iPO in dependentVars:
     2527                    depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
     2528                    if Ka2:
     2529                        depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
     2530            for i,name in enumerate(['omega','chi','phi']):
     2531                aname = pfx+'SH '+name
     2532                if aname in varylist:
     2533                    dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
     2534                    if Ka2:
     2535                        dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
     2536                elif aname in dependentVars:
     2537                    depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
     2538                    if Ka2:
     2539                        depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
     2540            for iSH in dFdODF:
     2541                if iSH in varylist:
     2542                    dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
     2543                    if Ka2:
     2544                        dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
     2545                elif iSH in dependentVars:
     2546                    depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
     2547                    if Ka2:
     2548                        depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
     2549            cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
     2550            for name,dpdA in cellDervNames:
     2551                if name in varylist:
     2552                    dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
     2553                    if Ka2:
     2554                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
     2555                elif name in dependentVars:
     2556                    depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
     2557                    if Ka2:
     2558                        depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
     2559            dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
     2560            for name in dDijDict:
     2561                if name in varylist:
     2562                    dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
     2563                    if Ka2:
     2564                        dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
     2565                elif name in dependentVars:
     2566                    depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
     2567                    if Ka2:
     2568                        depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
     2569            gamDict = GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict)
     2570            for name in gamDict:
     2571                if name in varylist:
     2572                    dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
     2573                    if Ka2:
     2574                        dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
     2575                elif name in dependentVars:
     2576                    depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
     2577                    if Ka2:
     2578                        depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
     2579                                               
     2580        elif 'T' in calcControls[hfx+'histType']:
     2581            print 'TOF Undefined at present'
     2582            raise Exception    #no TOF yet
     2583        #do atom derivatives -  for F,X & U so far             
     2584        corr = dervDict['int']/refl[9]
     2585        if Ka2:
     2586            corr2 = dervDict2['int']/refl[9]
     2587        for name in varylist+dependentVars:
     2588            try:
     2589                aname = name.split(pfx)[1][:2]
     2590                if aname not in ['Af','dA','AU']: continue # skip anything not an atom param
     2591            except IndexError:
     2592                continue
     2593            if name in varylist:
     2594                dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
     2595                if Ka2:
     2596                    dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
     2597            elif name in dependentVars:
     2598                depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
     2599                if Ka2:
     2600                    depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
     2601        # end loop over reflections
     2602    return depDerivDict,dMdv[:,iBegO:iFinO],iBegO,iFinO
     2603       
     2604def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):   
    23902605    # create a list of dependent variables and set up a dictionary to hold their derivatives
    23912606    dependentVars = G2mv.GetDependentVars()
     
    24132628        dx = x[1]-x[0]
    24142629        shl = max(parmDict[hfx+'SH/L'],0.002)
    2415         Ka2 = False
    24162630        if hfx+'Lam1' in parmDict.keys():
    24172631            wave = parmDict[hfx+'Lam1']
     
    24212635        else:
    24222636            wave = parmDict[hfx+'Lam']
     2637            Ka2 = False
     2638            lamRatio = None
     2639            kRatio = None
    24232640    else:
    24242641        print 'TOF Undefined at present'
     
    24362653        if 'Pawley' not in Phase['General']['Type']:
    24372654            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
    2438         for iref,refl in enumerate(refList):
    2439             if 'C' in calcControls[hfx+'histType']:        #CW powder
    2440                 h,k,l = refl[:3]
    2441                 dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
    2442                 if 'Pawley' in Phase['General']['Type']:
    2443                     try:
    2444                         refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
    2445                     except KeyError:
    2446 #                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
    2447                         continue
    2448                 Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
    2449                 iBeg = np.searchsorted(x,refl[5]-fmin)
    2450                 iFin = np.searchsorted(x,refl[5]+fmax)
    2451                 if not iBeg+iFin:       #peak below low limit - skip peak
     2655        else:
     2656            for iref,refl in enumerate(refList):
     2657                try:
     2658                    refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%refl[:3]])])
     2659                except KeyError:
     2660                    #print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%refl[:3]
    24522661                    continue
    2453                 elif not iBeg-iFin:     #peak above high limit - done
    2454                     break
    2455                 pos = refl[5]
    2456                 tanth = tand(pos/2.0)
    2457                 costh = cosd(pos/2.0)
    2458                 lenBF = iFin-iBeg
    2459                 dMdpk = np.zeros(shape=(6,lenBF))
    2460                 dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
    2461                 for i in range(1,5):
    2462                     dMdpk[i] += 100.*dx*refl[13]*refl[9]*dMdipk[i]
    2463                 dMdpk[0] += 100.*dx*refl[13]*refl[9]*dMdipk[0]
    2464                 dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
    2465                 if Ka2:
    2466                     pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
    2467                     kdelt = int((pos2-refl[5])/dx)               
    2468                     iBeg2 = min(lenX,iBeg+kdelt)
    2469                     iFin2 = min(lenX,iFin+kdelt)
    2470                     if iBeg2-iFin2:
    2471                         lenBF2 = iFin2-iBeg2
    2472                         dMdpk2 = np.zeros(shape=(6,lenBF2))
    2473                         dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
    2474                         for i in range(1,5):
    2475                             dMdpk2[i] = 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[i]
    2476                         dMdpk2[0] = 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[0]
    2477                         dMdpk2[5] = 100.*dx*refl[13]*dMdipk2[0]
    2478                         dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
    2479                 if 'Pawley' in Phase['General']['Type']:
    2480                     try:
    2481                         idx = varylist.index(pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)]))
    2482                         dMdv[idx][iBeg:iFin] = dervDict['int']/refl[9]
    2483                         if Ka2:
    2484                             dMdv[idx][iBeg2:iFin2] = dervDict2['int']/refl[9]
    2485                         # Assuming Pawley variables not in constraints
    2486                     except ValueError:
    2487                         pass
    2488                 dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
    2489                 names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
    2490                     hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
    2491                     hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
    2492                     hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
    2493                     hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
    2494                     hfx+'DisplaceY':[dpdY,'pos'],}
    2495                 for name in names:
    2496                     item = names[name]
    2497                     if name in varylist:
    2498                         dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
    2499                         if Ka2:
    2500                             dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
    2501                     elif name in dependentVars:
    2502                         if Ka2:
    2503                             depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
    2504                         depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
    2505                 for iPO in dIdPO:
    2506                     if iPO in varylist:
    2507                         dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
    2508                         if Ka2:
    2509                             dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
    2510                     elif iPO in dependentVars:
    2511                         depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
    2512                         if Ka2:
    2513                             depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
    2514                 for i,name in enumerate(['omega','chi','phi']):
    2515                     aname = pfx+'SH '+name
    2516                     if aname in varylist:
    2517                         dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
    2518                         if Ka2:
    2519                             dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
    2520                     elif aname in dependentVars:
    2521                         depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
    2522                         if Ka2:
    2523                             depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
    2524                 for iSH in dFdODF:
    2525                     if iSH in varylist:
    2526                         dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
    2527                         if Ka2:
    2528                             dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
    2529                     elif iSH in dependentVars:
    2530                         depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
    2531                         if Ka2:
    2532                             depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
    2533                 cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
    2534                 for name,dpdA in cellDervNames:
    2535                     if name in varylist:
    2536                         dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
    2537                         if Ka2:
    2538                             dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
    2539                     elif name in dependentVars:
    2540                         depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
    2541                         if Ka2:
    2542                             depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
    2543                 dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
    2544                 for name in dDijDict:
    2545                     if name in varylist:
    2546                         dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
    2547                         if Ka2:
    2548                             dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
    2549                     elif name in dependentVars:
    2550                         depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
    2551                         if Ka2:
    2552                             depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
    2553                 gamDict = GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict)
    2554                 for name in gamDict:
    2555                     if name in varylist:
    2556                         dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
    2557                         if Ka2:
    2558                             dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
    2559                     elif name in dependentVars:
    2560                         depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
    2561                         if Ka2:
    2562                             depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
    2563                                                
    2564             elif 'T' in calcControls[hfx+'histType']:
    2565                 print 'TOF Undefined at present'
    2566                 raise Exception    #no TOF yet
    2567             #do atom derivatives -  for F,X & U so far             
    2568             corr = dervDict['int']/refl[9]
    2569             if Ka2:
    2570                 corr2 = dervDict2['int']/refl[9]
    2571             for name in varylist+dependentVars:
    2572                 try:
    2573                     aname = name.split(pfx)[1][:2]
    2574                     if aname not in ['Af','dA','AU']: continue # skip anything not an atom param
    2575                 except IndexError:
    2576                     continue
    2577                 if name in varylist:
    2578                     dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
    2579                     if Ka2:
    2580                         dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
    2581                 elif name in dependentVars:
    2582                     depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
    2583                     if Ka2:
    2584                         depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
     2662        chunk = 200
     2663        if (calcControls['max Hprocess'] > 1 or
     2664            calcControls['max Rprocess'] == 1 or
     2665            mp.cpu_count() < 2 or len(refList) <  2 * chunk):
     2666            print 'single thread derivs'
     2667            TdepDerivDict, TdMdv, iBegO, iFinO = ComputeReflectionDerivative(
     2668                (0,refList,x,dependentVars,varylist,dFdvDict,
     2669                 Phase,calcControls,wave,G,g,GB,A,
     2670                 parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio))
     2671            dMdv[:,iBegO:iFinO] += TdMdv
     2672            for j in dependentVars:         
     2673                depDerivDict[j] += TdepDerivDict[j]
     2674        else:
     2675            print 'multiprocess derivs'
     2676            argList = []
     2677            for iref in range(0,len(refList),chunk):
     2678                argList.append(
     2679                    (iref,refList[iref:iref+chunk],x,dependentVars,varylist,dFdvDict,
     2680                     Phase,calcControls,wave,G,g,GB,A,
     2681                     parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio))
     2682            mpPool = mp.Pool()
     2683            for res in mpPool.imap_unordered(ComputeReflectionDerivative,argList):
     2684                if res is None: continue
     2685                TdepDerivDict, TdMdv, iBegO, iFinO = res
     2686                dMdv[:,iBegO:iFinO] += TdMdv
     2687                for j in dependentVars:
     2688                    depDerivDict[j] += TdepDerivDict[j]
     2689
    25852690    # now process derivatives in constraints
    25862691    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
     
    26752780                Histogram,parmdict,varylist,Phases,
    26762781                calcControls,pawleyLookup])
    2677     if MaxProcess > 1:
    2678         mpPool = mp.Pool(processes=MaxProcess)
    2679         #results = mpPool.map(ComputePowderHessian,argList)
    2680         #for Vec,Hess in results:
     2782    if MaxProcess > 1:
     2783        mpPool = mp.Pool()
    26812784        for Vec,Hess in mpPool.imap_unordered(ComputePowderHessian,argList):
    26822785            VecSum += Vec
    26832786            HessSum += Hess
     2787        mpPool.close()
     2788        mpPool.join()
    26842789    else:
    26852790        for arg in argList:
     
    27672872                )
    27682873    if MaxProcess > 1:
    2769         mpPool = mp.Pool(processes=MaxProcess)
     2874        mpPool = mp.Pool()
    27702875        for (xB,xF,ycSect,ybSect,RL,histkey
    27712876             ) in mpPool.imap_unordered(ComputePowderProfile,argList):
     
    27892894            Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
    27902895            M = np.concatenate((M,wdy))
     2896        mpPool.close()
     2897        mpPool.join()
    27912898    else:
    27922899        for arg in argList:
     
    33133420    print '\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]])   
    33143421
    3315            
    33163422def main():
    33173423    arg = sys.argv
Note: See TracChangeset for help on using the changeset viewer.