Changeset 2774


Ignore:
Timestamp:
Apr 7, 2017 4:05:51 PM (6 years ago)
Author:
vondreele
Message:

implement MC/SA for reflectometry
more image importer speed ups

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIplot.py

    r2772 r2774  
    28722872    def Draw():
    28732873        global xylim
    2874         print Page.Offset
    28752874        Plot.clear()
    28762875        Plot.set_title(Title)

  • trunk/GSASIIpwd.py

    r2772 r2774  
    19171917    print 'fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0])
    19181918   
     1919    class RandomDisplacementBounds(object):
     1920        """random displacement with bounds"""
     1921        def __init__(self, xmin, xmax, stepsize=0.5):
     1922            self.xmin = xmin
     1923            self.xmax = xmax
     1924            self.stepsize = stepsize
     1925   
     1926        def __call__(self, x):
     1927            """take a random step but ensure the new position is within the bounds"""
     1928            while True:
     1929                # this could be done in a much more clever way, but it will work for example purposes
     1930                steps = self.xmax-self.xmin
     1931                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
     1932                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
     1933                    break
     1934            return xnew
     1935   
    19191936    def GetModelParms():
    19201937        parmDict = {}
    … …  
    19221939        values = []
    19231940        bounds = []
    1924         parmDict['Res'] = data['Resolution'][0]/100.       #%-->decimal
     1941        parmDict['Res'] = data['Resolution'][0]/(100.*np.sqrt(ateln2))     #% FWHM-->decimal sig
    19251942        for parm in ['Scale','FltBack']:
    19261943            parmDict[parm] = data[parm][0]
    … …  
    19381955                if layer.get(parm,[0.,False])[1]:
    19391956                    varyList.append(cid+parm)
    1940                     values.append(layer[parm][0])
    1941                     bounds.append(Bounds[parm])
     1957                    value = layer[parm][0]
     1958                    values.append(value)
     1959                    if value:
     1960                        bound = [value*Bfac,value/Bfac]
     1961                    else:
     1962                        bound = [0.,10.]
     1963                    bounds.append(bound)
    19421964            if name not in ['vacuum','unit scatter']:
    19431965                parmDict[cid+'rho'] = Substances[name]['Scatt density']
    … …  
    20072029        Ic += (A**2+B**2)*Scale     
    20082030        return Ic
     2031       
     2032    def estimateT0(takestep):
     2033        Mmax = -1.e-10
     2034        Mmin = 1.e10
     2035        for i in range(100):
     2036            x0 = takestep(values)
     2037            M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList)
     2038            Mmin = min(M,Mmin)
     2039            MMax = max(M,Mmax)
     2040        return 1.5*(MMax-Mmin)
    20092041
    20102042    Q,Io,wt,Ic,Ib,Ifb = Profile[:6]
    … …  
    20142046    Qmax = Limits[1][1]
    20152047    wtFactor = ProfDict['wtFactor']
     2048    Bfac = data['Toler']
    20162049    Ibeg = np.searchsorted(Q,Qmin)
    20172050    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
    20182051    Ic[:] = 0
    2019     Bounds = {'Scale':[data['Scale'][0]*.85,data['Scale'][0]/.85],'FltBack':[None,None],
    2020               'DenMul':[None,None],'Thick':[1.,None],'Rough':[0.,None],'Mag SLD':[-10.,10.],'iDenMul':[None,None]}
     2052    Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6],
     2053              'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]}
    20212054    parmDict,varyList,values,bounds = GetModelParms()
    20222055    Msg = 'Failed to converge'
    … …  
    20302063            covM = result[1]
    20312064            newVals = result[0]
    2032         elif data['Minimizer'] == 'Global':
    2033             result = so.basinhopping(sumREFD,values,minimizer_kwargs={'method':'L-BFGS-B',
     2065        elif data['Minimizer'] == 'Basin Hopping':
     2066            xyrng = np.array(bounds).T
     2067            take_step = RandomDisplacementBounds(xyrng[0], xyrng[1])
     2068            T0 = estimateT0(take_step)
     2069            print ' Estimated temperature: %.3g'%(T0)
     2070            result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
     2071                interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
    20342072                'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList)})
    20352073            chisq = result.fun
    … …  
    20372075            newVals = result.x
    20382076            covM = []
     2077        elif data['Minimizer'] == 'MC/SA Anneal':
     2078            xyrng = np.array(bounds).T
     2079            result = G2mth.anneal(sumREFD, values, args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList),
     2080                schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
     2081                boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
     2082                ranRange=0.20,autoRan=False,dlg=None)
     2083            newVals = result[0]
     2084            parmDict.update(zip(varyList,newVals))
     2085            chisq = result[1]
     2086            ncalc = result[3]
     2087            covM = []
     2088            print ' MC/SA final temperature: %.4g'%(result[2])
    20392089        elif data['Minimizer'] == 'L-BFGS-B':
    20402090            result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds,   #ftol=Ftol,

  • trunk/GSASIIpwdGUI.py

    r2772 r2774  
    225225        {'Name':'vacuum','Rough':[0.,False],'Penetration':[0.,False],'DenMul':[1.0,False],}],   #bottom layer
    226226        'Scale':[1.0,False],'FltBack':[0.0,False],'Zero':'Top',                                 #globals
    227         'Minimizer':'LMLS','Resolution':[0.,'Const dq/q'],'Recomb':0.5,'Toler':0.001,           #minimizer controls
     227        'Minimizer':'LMLS','Resolution':[0.,'Const dq/q'],'Recomb':0.5,'Toler':0.5,           #minimizer controls
    228228        'DualFitFiles':['',],'DualFltBacks':[[0.0,False],],'DualScales':[[1.0,False],]}         #optional stuff for multidat fits?
    229229       
    … …  
    48534853        resol = wx.BoxSizer(wx.HORIZONTAL)
    48544854        resol.Add(wx.StaticText(G2frame.dataDisplay,label=' Instrument resolution (%'+GkDelta+'Q/Q): '),0,WACV)
    4855         resol.Add(G2G.ValidatedTxtCtrl(G2frame.dataDisplay,data['Resolution'],0,nDig=(10,2),typeHint=float),0,WACV)
     4855        resol.Add(G2G.ValidatedTxtCtrl(G2frame.dataDisplay,data['Resolution'],0,nDig=(10,2),min=0.,max=0.2),0,WACV)
    48564856        resol.Add(wx.StaticText(G2frame.dataDisplay,label=' Zero position location: '),0,WACV)
    48574857        poslist = ['Top','Bottom']
    … …  
    48634863        minimiz = wx.BoxSizer(wx.HORIZONTAL)
    48644864        minimiz.Add(wx.StaticText(G2frame.dataDisplay,label=' Minimizer: '),0,WACV)
    4865         minlist = ['LMLS','Global','L-BFGS-B',]
     4865        minlist = ['LMLS','Basin Hopping','MC/SA Anneal','L-BFGS-B',]
    48664866        minSel = wx.ComboBox(G2frame.dataDisplay,value=data['Minimizer'],choices=minlist,
    48674867            style=wx.CB_READONLY|wx.CB_DROPDOWN)
    48684868        minSel.Bind(wx.EVT_COMBOBOX, OnMinSel)
    48694869        minimiz.Add(minSel,0,WACV)
    4870         minimiz.Add(wx.StaticText(G2frame.dataDisplay,label=' Tolerance: '),0,WACV)
    4871         minimiz.Add(G2G.ValidatedTxtCtrl(G2frame.dataDisplay,data,'Toler',nDig=(10,1,'g'),typeHint=float),0,WACV)
     4870        minimiz.Add(wx.StaticText(G2frame.dataDisplay,label=' Bounds factor: '),0,WACV)
     4871        minimiz.Add(G2G.ValidatedTxtCtrl(G2frame.dataDisplay,data,'Toler',nDig=(10,2),max=0.99,min=0.1),0,WACV)
    48724872        weight = wx.CheckBox(G2frame.dataDisplay,label='Use 2% sig. weights')
    48734873        weight.SetValue(data.get('2% weight',False))

  • trunk/imports/G2img_1TIF.py

    r2773 r2774  
    208208                File.seek(8)
    209209                print 'Read GE-detector tiff file: ',filename
    210                 image = np.array(np.frombuffer(File.read(2*Npix),dtype=np.int16),dtype=np.int32)
     210                image = np.array(np.frombuffer(File.read(2*Npix),dtype=np.uint16),dtype=np.int32)
    211211#                image = np.fromfile(File,dtype=np.int16,count=2*Npix)[:Npix]
    212212#                image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
    … …  
    226226            File.seek(IFD[273][2][0])
    227227            print 'Read ImageJ tiff file: ',filename
    228             image = ar.array('H',File.read(2*Npix))
     228#            image = ar.array('H',File.read(2*Npix))
     229            image = File.read(2*Npix,dtype=np.uint16)
    229230            if '>' in byteOrd:
    230231                image.byteswap()
    231             image = np.array(np.asarray(image,dtype='H'),dtype=np.int32)           
     232            image = np.array(np.frombuffer(image),dtype=np.int32)
     233#            image = np.array(np.asarray(image,dtype='H'),dtype=np.int32)           
    232234    elif 262 in IFD and IFD[262][2][0] > 4:
    233235        tifType = 'DND'
    … …  
    235237        File.seek(512)
    236238        print 'Read DND SAX/WAX-detector tiff file: ',filename
    237         image = np.fromfile(File,dtype=np.int16,count=2*Npix)[:Npix]
     239        image = np.array(np.frombuffer(File.read(2*Npix),dtype=np.uint16),dtype=np.int32)
     240#        image = np.fromfile(File,dtype=np.int16,count=2*Npix)[:Npix]
    238241#        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
    239242    elif sizexy == [1536,1536]:
    … …  
    242245        File.seek(64)
    243246        print 'Read Gold tiff file:',filename
    244         image = np.fromfile(File,dtype=np.int16,count=2*Npix)[:Npix]
     247        image = np.array(np.frombuffer(File.read(2*Npix),dtype=np.uint16),dtype=np.int32)
     248#        image = np.fromfile(File,dtype=np.int16,count=2*Npix)[:Npix]
    245249#        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
    246250    elif sizexy == [2048,2048] or sizexy == [1024,1024] or sizexy == [3072,3072]:
    … …  
    256260#                    image = np.array(ar.array('f',File.read(4*Npix)),dtype=np.int32)
    257261                else:
    258                     image = np.fromfile(File,dtype=np.int,count=4*Npix)[:Npix]
     262                    image = np.array(np.frombuffer(File.read(4*Npix),dtype=np.int),dtype=np.int32)
     263#                    image = np.fromfile(File,dtype=np.int,count=4*Npix)[:Npix]
    259264#                    image = np.array(ar.array('I',File.read(4*Npix)),dtype=np.int32)
    260265            elif IFD[258][2][0] == 16:
    … …  
    263268                File.seek(8)
    264269                print 'Read MedOptics D1 tiff file: ',filename
    265                 image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
     270                image = np.array(np.frombuffer(File.read(2*Npix),dtype=np.uint16),dtype=np.int32)
     271#                image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
    266272                 
    267273        elif IFD[273][2][0] == 4096:
    … …  
    282288            File.seek(512)
    283289            print 'Read 11-ID-C tiff file: ',filename
    284             image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
     290            image = np.array(np.frombuffer(File.read(2*Npix),dtype=np.uint16),dtype=np.int32)
     291#            image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
    285292                   
    286293    elif sizexy == [4096,4096]:
    … …  
    297304                File.seek(8)
    298305                print 'Read PE 4Kx4K tiff file: ',filename
    299                 image = np.fromfile(File,dtype=np.uint,count=4*Npix)[:Npix]
    300                 if np.max(image) > 2**31:
    301                     image = np.array(image-2**30,dtype=np.int32)
    302                 else:
    303                     image = np.array(image,dtype=np.int32)
     306                image = np.array(np.frombuffer(File.read(4*Npix),dtype=np.float32)/2.**4,dtype=np.int32)
     307#                image = np.fromfile(File,dtype=np.uint,count=4*Npix)[:Npix]
     308#                if np.max(image) > 2**31:
     309#                    image = np.array(image-2**30,dtype=np.int32)
     310#                else:
     311#                    image = np.array(image,dtype=np.int32)
    304312#                arry = ar.array('I',File.read(4*Npix))
    305313#                image = np.array(arry)/2**16           
    … …  
    309317            File.seek(4096)
    310318            print 'Read Rayonix MX300HE tiff file: ',filename
    311             image = np.fromfile(File,dtype=np.int16,count=2*Npix)[:Npix]
     319            image = np.array(np.frombuffer(File.read(2*Npix),dtype=np.uint16),dtype=np.int32)
     320#            image = np.fromfile(File,dtype=np.int16,count=2*Npix)[:Npix]
    312321#            image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
    313322#    elif sizexy == [960,960]:
Note: See TracChangeset for help on using the changeset viewer.