Changeset 1252


Ignore:
Timestamp:
Mar 18, 2014 1:45:39 PM (10 years ago)
Author:
vondreele
Message:

unsuccessful implementation of IPG method (commented out) for size distribution
revisit later
plot SASD background on changes
document SASD shape form factors & volumes

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIplot.py

    r1250 r1252  
    890890                    D = xye[5]-Ymax*G2frame.delOffset
    891891                elif 'SASD' in plottype:
    892                     D = xye[4]-Ymax*G2frame.delOffset
     892                    D = xye[4]
    893893                    Plot.set_yscale("log",nonposy='mask')
    894894                    Plot.set_ylim(bottom=np.min(np.trim_zeros(Y))/2.,top=np.max(Y)*2.)
    … …  
    907907                        else:
    908908                            Plot.plot(X,Y,colors[N%6]+'+',picker=3.,clip_on=False)
     909                        Plot.plot(X,D,colors[(N+1)%6],picker=False)
    909910                        Plot.plot(X,Z,colors[(N+1)%6],picker=False)
    910911                elif G2frame.Weight and 'PWDR' in plottype:

  • trunk/GSASIIpwdGUI.py

    r1250 r1252  
    8686        'logBins':True,'Method':'MaxEnt','Distribution':[],
    8787        'Shape':['Spheroid',1.0],'MaxEnt':{'Niter':100,'Precision':0.01,'Sky':-3},
    88         'IPG':{'Niter':100,'Approach':0.8},'Reg':{},},           
     88        'IPG':{'Niter':100,'Approach':0.8,'Power':-1},'Reg':{},},           
    8989        'Unified':{'Levels':[],},           
    9090        'Particle':{'Levels':[],},
    … …  
    24752475    if isinstance(data['Size']['MaxEnt']['Sky'],float):
    24762476        data['Size']['MaxEnt']['Sky'] = -3
     2477    if 'Power' not in data['Size']['IPG']:
     2478        data['Size']['IPG']['Power'] = -1
    24772479    #end patches
    24782480   
    … …  
    24832485    def OnFitModel(event):
    24842486        print 'fit model for '+data['Current']
     2487        if not any(Sample['Contrast']):
     2488            G2frame.ErrorDialog('No contrast; your sample is a vacuum!',
     2489                'You need to define a scattering substance!\n'+    \
     2490                ' Do Substances and then Sample parameters')
     2491            return
    24852492        if data['Current'] == 'Size dist.':
    24862493            G2sasd.SizeDistribution(Profile,ProfDict,Limits,Substances,Sample,data)
    … …  
    25012508        Obj.SetValue(fmt%(value))
    25022509        data[itemKey][ind] = value
     2510        if itemKey == 'Back':
     2511            Profile[4][:] = value
     2512        G2plt.PlotPatterns(G2frame,plotType='SASD',newPlot=True)
    25032513       
    25042514    def OnCheckBox(event):
    … …  
    25872597        sizeSizer.Add((5,5),0)
    25882598        fitSizer = wx.BoxSizer(wx.HORIZONTAL)
    2589         methods = ['MaxEnt','IPG',]
     2599        methods = ['MaxEnt',]   #'IPG',]
    25902600        fitSizer.Add(wx.StaticText(G2frame.dataDisplay,label='Fitting method: '),0,WACV)
    25912601        method = wx.ComboBox(G2frame.dataDisplay,value=data['Size']['Method'],choices=methods,
    … …  
    26092619            floor.Bind(wx.EVT_COMBOBOX,OnIntVal)
    26102620            fitSizer.Add(floor,0,WACV)
     2621        elif 'IPG' in data['Size']['Method']:
     2622            fitSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Q power weight (-1 for sigma): '),0,WACV)
     2623            choices = ['-1','0','1','2','3','4']
     2624            power = wx.ComboBox(G2frame.dataDisplay,value=str(data['Size']['IPG']['Power']),choices=choices,
     2625                style=wx.CB_READONLY|wx.CB_DROPDOWN)
     2626            Indx[power.GetId()] = [data['Size']['IPG'],'Power',-2]
     2627            power.Bind(wx.EVT_COMBOBOX,OnIntVal)
     2628            fitSizer.Add(power,0,WACV)
    26112629        sizeSizer.Add(fitSizer,0)
    26122630
    … …  
    26682686    G2gd.HorizontalLine(mainSizer,G2frame.dataDisplay)
    26692687    if 'Size' in data['Current']:
     2688        if 'MaxEnt' in data['Size']['Method']:
     2689            Status.SetStatusText('Size distribution by Maximum entropy')
     2690        elif 'IPG' in data['Size']['Method']:
     2691            Status.SetStatusText('Size distribution by Interior-Point Gradient')
    26702692        mainSizer.Add(SizeSizer())       
    26712693    elif 'Particle' in data['Current']:
    … …  
    26812703    backVal.Bind(wx.EVT_KILL_FOCUS,OnValueChange)
    26822704    backSizer.Add(backVal,0,WACV)
    2683     backVar = wx.CheckBox(G2frame.dataDisplay,label='Apply?')
     2705    backVar = wx.CheckBox(G2frame.dataDisplay,label='Refine?')
    26842706    Indx[backVar.GetId()] = [data['Back'],1]
    26852707    backVar.SetValue(data['Back'][1])

  • trunk/GSASIIsasd.py

    r1250 r1252  
    6363def SphereFF(Q,R,args=()):
    6464    ''' Compute hard sphere form factor - can use numpy arrays
    65     param float:Q Q value array (usually in A-1)
    66     param float:R sphere radius (Usually in A - must match Q-1 units)
     65    param float Q: Q value array (usually in A-1)
     66    param float R: sphere radius (Usually in A - must match Q-1 units)
     67    param array args: ignored
    6768    returns float: form factors as array as needed
    6869    '''
    … …  
    7374    ''' Compute form factor of cylindrically symmetric ellipsoid (spheroid)
    7475    - can use numpy arrays for R & AR; will return corresponding numpy array
    75     param float:Q Q value array (usually in A-1)
     76    param float Q : Q value array (usually in A-1)
    7677    param float R: radius along 2 axes of spheroid
    77     param float AR: aspect ratio so 3rd axis = R*AR
     78    param array args: [float AR]: aspect ratio so 3rd axis = R*AR
    7879    returns float: form factors as array as needed
    7980    '''
    … …  
    8990def CylinderFF(Q,R,args):
    9091    ''' Compute form factor for cylinders - can use numpy arrays
    91     param float: Q Q value array (A-1)
    92     param float: R cylinder radius (A)
    93     param float: L cylinder length (A)
     92    param float Q: Q value array (A-1)
     93    param float R: cylinder radius (A)
     94    param array args: [float L]: cylinder length (A)
    9495    returns float: form factor
    9596    '''
    … …  
    107108def CylinderDFF(Q,L,args):
    108109    ''' Compute form factor for cylinders - can use numpy arrays
    109     param float: Q Q value array (A-1)
    110     param float: L cylinder half length (A)
    111     param float: R cylinder diameter (A)
     110    param float Q: Q value array (A-1)
     111    param float L: cylinder half length (A)
     112    param array args: [float R]: cylinder radius (A)
    112113    returns float: form factor
    113114    '''
    … …  
    117118def CylinderARFF(Q,R,args):
    118119    ''' Compute form factor for cylinders - can use numpy arrays
    119     param float: Q Q value array (A-1)
    120     param float: R cylinder radius (A)
    121     param float: AR cylinder aspect ratio = L/D = L/2R
     120    param float Q: Q value array (A-1)
     121    param float R: cylinder radius (A)
     122    param array args: [float AR]: cylinder aspect ratio = L/D = L/2R
    122123    returns float: form factor
    123124    '''
    … …  
    126127   
    127128def UniSphereFF(Q,R,args=0):
     129    ''' Compute form factor for unified sphere - can use numpy arrays
     130    param float Q: Q value array (A-1)
     131    param float R: cylinder radius (A)
     132    param array args: ignored
     133    returns float: form factor
     134    '''
    128135    Rg = np.sqrt(3./5.)*R
    129     B = 1.62/(Rg**4)    #are we missing *np.pi? 1.62 = 6*(3/5)**2/(4/3) sense?
     136    B = np.pi*1.62/(Rg**4)    #are we missing *np.pi? 1.62 = 6*(3/5)**2/(4/3) sense?
    130137    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg/np.sqrt(6)))**3
    131138    return np.sqrt(np.exp((-Q[:,np.newaxis]**2*Rg**2)/3.)+(B/QstV**4))
    132139   
    133140def UniRodFF(Q,R,args):
     141    ''' Compute form factor for unified rod - can use numpy arrays
     142    param float Q: Q value array (A-1)
     143    param float R: cylinder radius (A)
     144    param array args: [float R]: cylinder radius (A)
     145    returns float: form factor
     146    '''
    134147    L = args[0]
    135148    Rg2 = np.sqrt(R**2/2+L**2/12)
    … …  
    145158   
    146159def UniRodARFF(Q,R,args):
     160    ''' Compute form factor for unified rod of fixed aspect ratio - can use numpy arrays
     161    param float Q: Q value array (A-1)
     162    param float R: cylinder radius (A)
     163    param array args: [float AR]: cylinder aspect ratio = L/D = L/2R
     164    returns float: form factor
     165    '''
    147166    AR = args[0]
    148167    return UniRodFF(Q,R,[AR*R,])
    149168   
    150169def UniDiskFF(Q,R,args):
     170    ''' Compute form factor for unified disk - can use numpy arrays
     171    param float Q: Q value array (A-1)
     172    param float R: cylinder radius (A)
     173    param array args: [float T]: disk thickness (A)
     174    returns float: form factor
     175    '''
    151176    T = args[0]
    152177    Rg2 = np.sqrt(R**2/2.+T**2/12.)
    … …  
    163188   
    164189def UniTubeFF(Q,R,args):
     190    ''' Compute form factor for unified disk - can use numpy arrays
     191    param float Q: Q value array (A-1)
     192    param float R: cylinder radius (A)
     193    param array args: [float L,T]: tube length & wall thickness(A)
     194    returns float: form factor
     195    '''
    165196    L,T = args[:2]
    166197    Ri = R-T
    … …  
    186217    ''' Compute volume of sphere
    187218    - numpy array friendly
    188     param float:R sphere radius
     219    param float R: sphere radius
     220    param array args: ignored
    189221    returns float: volume
    190222    '''
    … …  
    195227    - numpy array friendly
    196228    param float R: radius along 2 axes of spheroid
    197     param float AR: aspect ratio so radius of 3rd axis = R*AR
     229    param array args: [float AR]: aspect ratio so radius of 3rd axis = R*AR
    198230    returns float: volume
    199231    '''
    … …  
    204236    ''' Compute cylinder volume for radius & length
    205237    - numpy array friendly
    206     param float: R diameter (A)
    207     param float: L length (A)
     238    param float R: diameter (A)
     239    param array args: [float L]: length (A)
    208240    returns float:volume (A^3)
    209241    '''
    … …  
    215247    - numpy array friendly
    216248    param float: L half length (A)
    217     param float: D diameter (A)
     249    param array args: [float D]: diameter (A)
    218250    returns float:volume (A^3)
    219251    '''
    … …  
    225257    - numpy array friendly
    226258    param float: R radius (A)
    227     param float: AR=L/D=L/2R aspect ratio
     259    param array args: [float AR]: =L/D=L/2R aspect ratio
    228260    returns float:volume
    229261    '''
    … …  
    234266    ''' Compute volume of sphere
    235267    - numpy array friendly
    236     param float:R sphere radius
     268    param float R: sphere radius
     269    param array args: ignored
    237270    returns float: volume
    238271    '''
    … …  
    242275    ''' Compute cylinder volume for radius & length
    243276    - numpy array friendly
    244     param float: R diameter (A)
    245     param float: L length (A)
     277    param float R: diameter (A)
     278    param array args: [float L]: length (A)
    246279    returns float:volume (A^3)
    247280    '''
    … …  
    250283   
    251284def UniRodARVol(R,args):
     285    ''' Compute rod volume for radius & aspect ratio
     286    - numpy array friendly
     287    param float R: diameter (A)
     288    param array args: [float AR]: =L/D=L/2R aspect ratio
     289    returns float:volume (A^3)
     290    '''
    252291    AR = args[0]
    253292    return CylinderARVol(R,[AR,])
    254293   
    255294def UniDiskVol(R,args):
     295    ''' Compute disk volume for radius & thickness
     296    - numpy array friendly
     297    param float R: diameter (A)
     298    param array args: [float T]: thickness
     299    returns float:volume (A^3)
     300    '''
    256301    T = args[0]
    257302    return CylinderVol(R,[T,])
    … …  
    260305    ''' Compute tube volume for radius, length & wall thickness
    261306    - numpy array friendly
    262     param float: R diameter (A)
    263     param float: L length (A)
    264     param float: T tube wall thickness (A)
     307    param float R: diameter (A)
     308    param array args: [float L,T]: tube length & wall thickness(A)
    265309    returns float: volume (A^3) of tube wall
    266310    '''
    … …  
    317361    pass
    318362
    319 def MaxEnt_SB(datum, sigma, base, IterMax, G, image_to_data=None, data_to_image=None, report=False):
     363def MaxEnt_SB(datum, sigma, G, base, IterMax, image_to_data=None, data_to_image=None, report=False):
    320364    '''
    321365    do the complete Maximum Entropy algorithm of Skilling and Bryan
    … …  
    323367    :param float datum[]:
    324368    :param float sigma[]:
     369    :param float[][] G: transformation matrix
    325370    :param float base[]:
    326371    :param int IterMax:
    327     :param float[][] G: transformation matrix
    328372    :param obj image_to_data: opus function (defaults to opus)
    329373    :param obj data_to_image: tropus function (defaults to tropus)
    … …  
    617661    return chisq,f,image_to_data(f, G)       # no solution after IterMax iterations
    618662
     663   
     664###############################################################################
     665#### IPG/TNNLS Routines
     666###############################################################################
     667
     668def IPG(datum,sigma,G,Bins,Dbins,IterMax,Qvec=[],approach=0.8,Power=-1,report=False):
     669    ''' An implementation of the Interior-Point Gradient method of
     670    Michael Merritt & Yin Zhang, Technical Report TR04-08, Dept. of Comp. and
     671    Appl. Math., Rice Univ., Houston, Texas 77005, U.S.A. found on the web at
     672    http://www.caam.rice.edu/caam/trs/2004/TR04-08.pdf
     673    Problem addressed: Total Non-Negative Least Squares (TNNLS)
     674    :param float datum[]:
     675    :param float sigma[]:
     676    :param float[][] G: transformation matrix
     677    :param int IterMax:
     678    :param float Qvec: data positions for Power = 0-4
     679    :param float approach: 0.8 default fitting parameter
     680    :param int Power: 0-4 for Q^Power weighting, -1 to use input sigma
     681   
     682    '''
     683    if Power < 0:
     684        GmatE = G/sigma[:np.newaxis]
     685        IntE = datum/sigma
     686        pwr = 0
     687        QvecP = np.ones_like(datum)
     688    else:
     689        GmatE = G[:]
     690        IntE = datum[:]
     691        pwr = Power
     692        QvecP = Qvec**pwr
     693    Amat = GmatE*QvecP[:np.newaxis]
     694    AAmat = np.inner(Amat,Amat)
     695    Bvec = datum*QvecP
     696    Xw = np.ones_like(Bins)*1.e-6
     697    calc = np.dot(G.T,Xw)
     698    nIter = 0
     699    err = 10.
     700    while (nIter<IterMax) and (err > 1.):
     701        #Step 1 in M&Z paper:
     702        Qk = np.dot(AAmat,Xw)-np.dot(Amat,Bvec)
     703        Dk = Xw/np.dot(AAmat,Xw)
     704        Pk = -Dk*Qk
     705        #Step 2 in M&Z paper:
     706        alpSt = -np.dot(Pk,Qk)/np.dot(Pk,np.dot(AAmat,Pk))
     707        alpWv = -Xw/Pk
     708        AkSt = [np.where(Pk[i]<0,np.min((approach*alpWv[i],alpSt)),Pk[i]) for i in range(Pk.shape[0])]
     709        #Step 3 in M&Z paper:
     710        shift = AkSt*Pk
     711        print np.sum(shift**2)
     712        Xw += shift
     713        #done IPG; now results
     714        nIter += 1
     715        calc = np.dot(G.T,Xw)
     716        chisq = np.sum(((datum-calc)/sigma)**2)
     717        err = chisq/len(datum)
     718        if report:
     719            print ' Iteration: %d, chisq: %.3g'%(nIter,chisq)
     720    return chisq,Xw,calc
     721
     722###############################################################################
     723#### SASD Utilities
     724###############################################################################
     725
     726def SetScale(Data,refData):
     727    Profile,Limits,Sample = Data
     728    refProfile,refLimits,refSample = refData
     729    x,y = Profile[:2]
     730    rx,ry = refProfile[:2]
     731    Beg = np.max([rx[0],x[0],Limits[1][0],refLimits[1][0]])
     732    Fin = np.min([rx[-1],x[-1],Limits[1][1],refLimits[1][1]])
     733    iBeg = np.searchsorted(x,Beg)
     734    iFin = np.searchsorted(x,Fin)
     735    sum = np.sum(y[iBeg:iFin])
     736    refsum = np.sum(np.interp(x[iBeg:iFin],rx,ry,0,0))
     737    Sample['Scale'][0] = refSample['Scale'][0]*refsum/sum
     738   
     739###############################################################################
     740#### Size distribution
     741###############################################################################
     742
     743def SizeDistribution(Profile,ProfDict,Limits,Substances,Sample,data):
     744    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
     745        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
     746        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
     747        'Unified disk':[UniDiskFF,UniDiskVol]}
     748    Shape = data['Size']['Shape'][0]
     749    Parms = data['Size']['Shape'][1:]
     750    if data['Size']['logBins']:
     751        Bins = np.logspace(np.log10(data['Size']['MinDiam']),np.log10(data['Size']['MaxDiam']),
     752            data['Size']['Nbins']+1,True)/2.        #make radii
     753    else:
     754        Bins = np.linspace(data['Size']['MinDiam'],data['Size']['MaxDiam'],
     755            data['Size']['Nbins']+1,True)/2.        #make radii
     756    Dbins = np.diff(Bins)
     757    Bins = Bins[:-1]+Dbins/2.
     758    Contrast = Sample['Contrast'][1]
     759    Scale = Sample['Scale'][0]
     760    Sky = 10**data['Size']['MaxEnt']['Sky']
     761    BinsBack = np.ones_like(Bins)*Sky*Scale/Contrast #How about *Scale/Contrast?
     762    Back = data['Back']
     763    Q,Io,wt,Ic,Ib = Profile
     764    Qmin = Limits[1][0]
     765    Qmax = Limits[1][1]
     766    wtFactor = ProfDict['wtFactor']
     767    Ibeg = np.searchsorted(Q,Qmin)
     768    Ifin = np.searchsorted(Q,Qmax)
     769    BinMag = np.zeros_like(Bins)
     770    Ic[:] = 0.
     771    Gmat = G_matrix(Q[Ibeg:Ifin],Bins,Contrast,shapes[Shape][0],shapes[Shape][1],args=Parms)
     772    if 'MaxEnt' == data['Size']['Method']:
     773        chisq,BinMag,Ic[Ibeg:Ifin] = MaxEnt_SB(Scale*Io[Ibeg:Ifin]-Back[0],
     774            Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),Gmat,BinsBack,
     775            data['Size']['MaxEnt']['Niter'],report=True)
     776    elif 'IPG' == data['Size']['Method']:
     777        chisq,BinMag,Ic[Ibeg:Ifin] = IPG(Scale*Io[Ibeg:Ifin]-Back[0],Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),
     778            Gmat,Bins,Dbins,data['Size']['IPG']['Niter'],Q[Ibeg:Ifin],approach=0.8,
     779            Power=data['Size']['IPG']['Power'],report=True)
     780    Ib[:] = Back[0]
     781    Ic[Ibeg:Ifin] += Back[0]
     782    print ' Final chi^2: %.3f'%(chisq)
     783    Vols = shapes[Shape][1](Bins,Parms)
     784    data['Size']['Distribution'] = [Bins,Dbins,BinMag/(2.*Dbins)]
     785       
    619786   
    620787################################################################################
    … …  
    685852if __name__ == '__main__':
    686853    tests()
    687    
    688 ###############################################################################
    689 #### SASD Utilities
    690 ###############################################################################
    691 
    692 def SetScale(Data,refData):
    693     Profile,Limits,Sample = Data
    694     refProfile,refLimits,refSample = refData
    695     x,y = Profile[:2]
    696     rx,ry = refProfile[:2]
    697     Beg = np.max([rx[0],x[0],Limits[1][0],refLimits[1][0]])
    698     Fin = np.min([rx[-1],x[-1],Limits[1][1],refLimits[1][1]])
    699     iBeg = np.searchsorted(x,Beg)
    700     iFin = np.searchsorted(x,Fin)
    701     sum = np.sum(y[iBeg:iFin])
    702     refsum = np.sum(np.interp(x[iBeg:iFin],rx,ry,0,0))
    703     Sample['Scale'][0] = refSample['Scale'][0]*refsum/sum
    704    
    705 ###############################################################################
    706 #### Size distribution
    707 ###############################################################################
    708 
    709 def SizeDistribution(Profile,ProfDict,Limits,Substances,Sample,data):
    710     shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
    711         'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
    712         'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
    713         'Unified disk':[UniDiskFF,UniDiskVol]}
    714     Shape = data['Size']['Shape'][0]
    715     Parms = data['Size']['Shape'][1:]
    716     if data['Size']['logBins']:
    717         Bins = np.logspace(np.log10(data['Size']['MinDiam']),np.log10(data['Size']['MaxDiam']),
    718             data['Size']['Nbins']+1,True)/2.        #make radii
    719     else:
    720         Bins = np.linspace(data['Size']['MinDiam'],data['Size']['MaxDiam'],
    721             data['Size']['Nbins']+1,True)/2.        #make radii
    722     Dbins = np.diff(Bins)
    723     Bins = Bins[:-1]+Dbins/2.
    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?
    728     Back = data['Back']
    729     Q,Io,wt,Ic,Ib = Profile[:5]
    730     Qmin = Limits[1][0]
    731     Qmax = Limits[1][1]
    732     wtFactor = ProfDict['wtFactor']
    733     Ibeg = np.searchsorted(Q,Qmin)
    734     Ifin = np.searchsorted(Q,Qmax)
    735     Gmat = G_matrix(Q[Ibeg:Ifin],Bins,Contrast,shapes[Shape][0],shapes[Shape][1],args=Parms)
    736     chisq,BinMag,Ic[Ibeg:Ifin] = MaxEnt_SB(Scale*Io[Ibeg:Ifin]-Back[0],
    737         Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),BinsBack,
    738         data['Size']['MaxEnt']['Niter'],Gmat,report=True)
    739     if Back[1]:
    740         Ib = Back[0]
    741         Ic[Ibeg:Ifin] += Back[0]
    742     print ' Final chi^2: %.3f'%(chisq)
    743     Vols = shapes[Shape][1](Bins,Parms)
    744     data['Size']['Distribution'] = [Bins,Dbins,BinMag/(2.*Dbins)]
    745        
    746    
     854
Note: See TracChangeset for help on using the changeset viewer.