Changeset 1249


Ignore:
Timestamp:
Mar 14, 2014 11:18:21 AM (10 years ago)
Author:
vondreele
Message:

include the 'sky' parameter in MaxEnt?; scale it by Scale & Contrast - makes it more universal.
Reorder the opus & tropus routines back to original definitions as per Pete J.

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIpwdGUI.py

    r1248 r1249  
    8585    return {'Back':[0.0,False],'Size':{'MinDiam':50,'MaxDiam':10000,'Nbins':100,
    8686        'logBins':True,'Method':'MaxEnt','Distribution':[],
    87         'Shape':['Spheroid',1.0],'MaxEnt':{'Niter':100,'Precision':0.01,'Sky':1e-6},
     87        'Shape':['Spheroid',1.0],'MaxEnt':{'Niter':100,'Precision':0.01,'Sky':-6},
    8888        'IPG':{'Niter':100,'Approach':0.8},'Reg':{},},           
    8989        'Unified':{'Levels':[],},           
    … …  
    24732473        data['Size']['MaxDiam'] = 10000.
    24742474        del data['Size']['MinMaxDiam']
     2475    if isinstance(data['Size']['MaxEnt']['Sky'],float):
     2476        data['Size']['MaxEnt']['Sky'] = -3
    24752477    #end patches
    24762478   
    … …  
    25072509    def OnIntVal(event):
    25082510        Obj = event.GetEventObject()
    2509         item,ind = Indx[Obj.GetId()]
    2510         item[ind] = int(Obj.GetValue())
    2511        
     2511        item,ind,minVal = Indx[Obj.GetId()]
     2512        try:
     2513            value = int(Obj.GetValue())
     2514            if value <= minVal:
     2515                raise ValueError
     2516        except ValueError:
     2517            value = item[ind]
     2518            Obj.SetValue(str(value))
     2519        item[ind] = value
     2520
    25122521    def SizeSizer():
    25132522       
    … …  
    25352544        nbins = wx.ComboBox(G2frame.dataDisplay,value=str(data['Size']['Nbins']),choices=bins,
    25362545            style=wx.CB_READONLY|wx.CB_DROPDOWN)
    2537         Indx[nbins.GetId()] = [data['Size'],'Nbins']
     2546        Indx[nbins.GetId()] = [data['Size'],'Nbins',0]
    25382547        nbins.Bind(wx.EVT_COMBOBOX,OnIntVal)       
    25392548        binSizer.Add(nbins,0,WACV)
    … …  
    25432552            style=wx.CB_READONLY|wx.CB_DROPDOWN)
    25442553        mindiam.Bind(wx.EVT_COMBOBOX,OnIntVal)       
    2545         Indx[mindiam.GetId()] = [data['Size'],'MinDiam']
     2554        Indx[mindiam.GetId()] = [data['Size'],'MinDiam',0]
    25462555        binSizer.Add(mindiam,0,WACV)
    25472556        binSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Max diam.: '),0,WACV)
    25482557        maxDias = [str(1000*(i+1)) for i in range(10)]
    25492558        maxdiam = wx.ComboBox(G2frame.dataDisplay,value=str(data['Size']['MaxDiam']),choices=maxDias,
    2550             style=wx.CB_READONLY|wx.CB_DROPDOWN)
     2559            style=wx.CB_DROPDOWN)
    25512560        maxdiam.Bind(wx.EVT_COMBOBOX,OnIntVal)       
    2552         Indx[maxdiam.GetId()] = [data['Size'],'MaxDiam']
     2561        Indx[maxdiam.GetId()] = [data['Size'],'MaxDiam',0]
    25532562        binSizer.Add(maxdiam,0,WACV)
    25542563        logbins = wx.CheckBox(G2frame.dataDisplay,label='Log bins?')
    … …  
    25892598        iter = wx.ComboBox(G2frame.dataDisplay,value=str(data['Size'][Method]['Niter']),choices=iters,
    25902599            style=wx.CB_READONLY|wx.CB_DROPDOWN)
    2591         Indx[iter.GetId()] = [data['Size'][Method],'Niter']
     2600        Indx[iter.GetId()] = [data['Size'][Method],'Niter',0]
    25922601        iter.Bind(wx.EVT_COMBOBOX,OnIntVal)
    25932602        fitSizer.Add(iter,0,WACV)
    2594        
     2603        if 'MaxEnt' in data['Size']['Method']:
     2604            fitSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Log floor factor: '),0,WACV)
     2605            floors = [str(-i) for i in range(9)]
     2606            floor = wx.ComboBox(G2frame.dataDisplay,value=str(data['Size']['MaxEnt']['Sky']),choices=floors,
     2607                style=wx.CB_READONLY|wx.CB_DROPDOWN)
     2608            Indx[floor.GetId()] = [data['Size']['MaxEnt'],'Sky',-10]
     2609            floor.Bind(wx.EVT_COMBOBOX,OnIntVal)
     2610            fitSizer.Add(floor,0,WACV)
    25952611        sizeSizer.Add(fitSizer,0)
    25962612
    … …  
    26602676    backSizer = wx.BoxSizer(wx.HORIZONTAL)
    26612677    backSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Background:'),0,WACV)
    2662     backVal = wx.TextCtrl(G2frame.dataDisplay,value='%.3f'%(data['Back'][0]),style=wx.TE_PROCESS_ENTER)
    2663     Indx[backVal.GetId()] = ['Back',0,'%.3f']
     2678    backVal = wx.TextCtrl(G2frame.dataDisplay,value='%.3g'%(data['Back'][0]),style=wx.TE_PROCESS_ENTER)
     2679    Indx[backVal.GetId()] = ['Back',0,'%.3g']
    26642680    backVal.Bind(wx.EVT_TEXT_ENTER,OnValueChange)       
    26652681    backVal.Bind(wx.EVT_KILL_FOCUS,OnValueChange)

  • trunk/GSASIIsasd.py

    r1248 r1249  
    1313# $Id: GSASIIsasd.py 1186 2014-01-09 17:09:53Z vondreele $
    1414########### SVN repository information ###################
     15import os
    1516import sys
    1617import math
    … …  
    283284    FF = FFfxn(q,r,args)
    284285    Vol = Volfxn(r,args)
    285     return 1.e-4*(contrast*Vol*FF**2)     #10^-20 vs 10^-24
     286    return 1.e-4*(contrast*Vol*FF**2).T     #10^-20 vs 10^-24
    286287   
    287288'''
    … …  
    312313'''
    313314
    314 import os
    315 import sys
    316 import math
    317 import numpy
    318 
    319315class MaxEntException(Exception):
    320316    '''Any exception from this module'''
    … …  
    345341    MOVE_PASSES       = 0.001                   # convergence test in routine: move
    346342
    347     def opus (data, G):
    348         '''
    349         opus: transform data-space -> solution-space:  [G] * data
     343    def tropus (data, G):
     344        '''
     345        tropus: transform data-space -> solution-space:  [G] * data
    350346       
    351347        default definition, caller can use this definition or provide an alternative
    … …  
    357353        return G.dot(data)
    358354
    359     def tropus (image, G):
    360         '''
    361         tropus: transform solution-space -> data-space:  [G]^tr * image
     355    def opus (image, G):
     356        '''
     357        opus: transform solution-space -> data-space:  [G]^tr * image
    362358       
    363359        default definition, caller can use this definition or provide an alternative
    … …  
    367363        :returns float[M]: calculated data, ndarray of shape (M)
    368364        '''
    369         return G.transpose().dot(image)
     365        return np.dot(G.T,image)    #G.transpose().dot(image)
    370366
    371367    def Dist(s2, beta):
    … …  
    402398        '''
    403399        n = b.shape[0]
    404         fl = numpy.ndarray((n, n))*0
    405         bl = numpy.ndarray((n))*0
     400        fl = np.zeros((n,n))
     401        bl = np.zeros_like(b)
    406402       
    407403        #print_arr("ChoSol: a", a)
    … …  
    443439   
    444440        # last, compute beta from bl and fl
    445         beta = numpy.ndarray((n))
     441        beta = np.empty((n))
    446442        beta[-1] = bl[-1] / fl[-1][-1]
    447443        for i in (1, 0):
    … …  
    511507    # to enable parts of them to be passed as
    512508    # as vectors to "image_to_data" and "data_to_image".
    513     xi      = 0*numpy.ndarray((SEARCH_DIRECTIONS, n))
    514     eta     = 0*numpy.ndarray((SEARCH_DIRECTIONS, npt))
    515     beta    = 0*numpy.ndarray((SEARCH_DIRECTIONS))
    516     # s1      = 0*numpy.ndarray((SEARCH_DIRECTIONS))
    517     # c1      = 0*numpy.ndarray((SEARCH_DIRECTIONS))
    518     s2      = 0*numpy.ndarray((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
    519     c2      = 0*numpy.ndarray((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
     509    xi      = np.zeros((SEARCH_DIRECTIONS, n))
     510    eta     = np.zeros((SEARCH_DIRECTIONS, npt))
     511    beta    = np.zeros((SEARCH_DIRECTIONS))
     512    # s1      = np.zeros((SEARCH_DIRECTIONS))
     513    # c1      = np.zeros((SEARCH_DIRECTIONS))
     514    s2      = np.zeros((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
     515    c2      = np.zeros((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
    520516
    521517    # TODO: replace blank (scalar) with base (vector)
    … …  
    533529
    534530        cgrad = data_to_image (ox, G)              # cgrad[i] = del(C)/del(f[i]), SB eq. 8
    535         sgrad = -numpy.log(f/base) / (blank*math.exp (1.0))  # sgrad[i] = del(S)/del(f[i])
     531        sgrad = -np.log(f/base) / (blank*math.exp (1.0))  # sgrad[i] = del(S)/del(f[i])
    536532        snorm = math.sqrt(sum(f * sgrad*sgrad))    # entropy term, SB eq. 22
    537533        cnorm = math.sqrt(sum(f * cgrad*cgrad))    # ChiSqr term, SB eq. 22
    … …  
    599595
    600596        # calculate the normalized entropy
    601         S = sum((f/fSum) * numpy.log(f/fSum))      # normalized entropy, S&B eq. 1
     597        S = sum((f/fSum) * np.log(f/fSum))      # normalized entropy, S&B eq. 1
    602598        z = (datum - image_to_data (f, G)) / sigma  # standardized residuals
    603599        chisq = sum(z*z)                            # report this ChiSq
    … …  
    654650        M = len(buf)
    655651        buf = zip(*buf)         # transpose rows and columns
    656         q  = numpy.array(buf[0], dtype=numpy.float64)
    657         I  = numpy.array(buf[1], dtype=numpy.float64)
    658         dI = numpy.array(buf[2], dtype=numpy.float64)
     652        q  = np.array(buf[0], dtype=np.float64)
     653        I  = np.array(buf[1], dtype=np.float64)
     654        dI = np.array(buf[2], dtype=np.float64)
    659655        return q, I, dI
    660656    print "MaxEnt_SB: "
    … …  
    667663    errFac = 1.05
    668664   
    669     r    = numpy.logspace(math.log10(dMin), math.log10(dMax), nRadii)/2
     665    r    = np.logspace(math.log10(dMin), math.log10(dMax), nRadii)/2
    670666    dr   = r * (r[1]/r[0] - 1)          # step size
    671     f_dr = numpy.ndarray((nRadii)) * 0  # volume fraction histogram
    672     b    = numpy.ndarray((nRadii)) * 0 + defaultDistLevel  # MaxEnt "sky background"
     667    f_dr = np.ndarray((nRadii)) * 0  # volume fraction histogram
     668    b    = np.ndarray((nRadii)) * 0 + defaultDistLevel  # MaxEnt "sky background"
    673669   
    674670    qVec, I, dI = readTextData(test_data_file)
    … …  
    726722    Dbins = np.diff(Bins)
    727723    Bins = Bins[:-1]+Dbins/2.
    728     BinsBack = np.ones_like(Bins)*1.e-6
     724    Contrast = Sample['Contrast'][1]
     725    Scale = Sample['Scale'][0]
     726    Sky = 10**data['Size']['MaxEnt']['Sky']
     727    BinsBack = np.ones_like(Bins)*Sky*Scale/Contrast #How about *Scale/Contrast?
    729728    Back = data['Back']
    730729    Q,Io,wt,Ic,Ib = Profile[:5]
    731730    Qmin = Limits[1][0]
    732731    Qmax = Limits[1][1]
    733     Contrast = Sample['Contrast'][1]
    734732    wtFactor = ProfDict['wtFactor']
    735733    Ibeg = np.searchsorted(Q,Qmin)
    … …  
    739737        Ic[Ibeg:Ifin] = Back[0]
    740738    Gmat = G_matrix(Q[Ibeg:Ifin],Bins,Contrast,shapes[Shape][0],shapes[Shape][1],args=Parms)
    741     chisq,BinMag,Ic[Ibeg:Ifin] = MaxEnt_SB(Sample['Scale'][0]*Io[Ibeg:Ifin]-Back[0],
    742         Sample['Scale'][0]/np.sqrt(wtFactor*wt[Ibeg:Ifin]),BinsBack,
     739    chisq,BinMag,Ic[Ibeg:Ifin] = MaxEnt_SB(Scale*Io[Ibeg:Ifin]-Back[0],
     740        Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),BinsBack,
    743741        data['Size']['MaxEnt']['Niter'],Gmat,report=True)
    744742    print ' Final chi^2: %.3f'%(chisq)
Note: See TracChangeset for help on using the changeset viewer.