Changeset 4094


Ignore:
Timestamp:
Aug 15, 2019 2:54:58 PM (2 years ago)
Author:
vondreele
Message:

SHAPES now complete, but slow.

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/G2shapes.py

    r4084 r4094  
    2525import sys
    2626import os
     27import copy
    2728import random
    2829import time
     
    647648   
    648649        angstrom_scale = 1.0
    649         Bins,Dbins,BinMag = data['Size']['Distribution']
     650        Bins,Dbins,BinMag = data['Pair']['Distribution']
    650651       
    651652        aList_r += list(Bins)
     
    18301831        # Write input and model P(r)
    18311832#        pr_writer(aList_pr,aList_r,aList_pr_model,file_pr)
    1832         PRcalc.append([aList_r,aList_pr,aList_pr_model,delta_hist_sum])
     1833        PRcalc.append([aList_r,aList_pr,copy.copy(aList_pr_model),delta_hist_sum])
    18331834   
    18341835        # Calculate comparison versus intensities
  • trunk/GSASIIplot.py

    r4092 r4094  
    91519151    Gr = np.array([0,255,0])
    91529152    Bl = np.array([0,0,255])
    9153     uBox = np.array([[0,0,0],[1,0,0],[0,1,0],[0,0,1]])
     9153    uBox = np.array([[0,0,0],[50,0,0],[0,50,0],[0,0,50]])
    91549154    uEdges = np.array([[uBox[0],uBox[1]],[uBox[0],uBox[2]],[uBox[0],uBox[3]]])
    91559155    uColors = [Rd,Gr,Bl]
     
    92889288        GL.glLoadIdentity()
    92899289        GL.glViewport(0,0,VS[0],VS[1])
    9290         GLU.gluPerspective(20.,aspect,1.,500.)
     9290        GLU.gluPerspective(50.,aspect,1.,500.)
    92919291        GLU.gluLookAt(0,0,cPos,0,0,0,0,1,0)
    92929292        SetLights()           
  • trunk/GSASIIpwdGUI.py

    r4092 r4094  
    356356                'MaxEnt':{'Niter':100,'Precision':0.01,'Sky':-3},
    357357                'IPG':{'Niter':100,'Approach':0.8,'Power':-1},'Reg':{},},
    358         'Pair':{'Method':'Regularization','MaxRadius':100.,'NBins':100,'Errors':'User',
     358        'Pair':{'Method':'Moore','MaxRadius':100.,'NBins':100,'Errors':'User',
    359359                'Percent error':2.5,'Background':[0,False],'Distribution':[],
    360360                'Moore':20,'Dist G':100.,'Result':[],},           
     
    52305230        data['BackFile'] = ''
    52315231    if 'Pair' not in data:
    5232         data['Pair'] = {'Method':'Regularization','MaxRadius':100.,'NBins':100,'Errors':'User','Result':[],
     5232        data['Pair'] = {'Method':'Moore','MaxRadius':100.,'NBins':100,'Errors':'User','Result':[],
    52335233            'Percent error':2.5,'Background':[0,False],'Distribution':[],'Moore':[20,False],'Dist G':100.,} 
    52345234    if 'Shapes' not in data:
     
    54225422            SaveState()
    54235423            G2sasd.PairDistFxn(Profile,ProfDict,Limits,Sample,data)
     5424            RefreshPlots(True)
    54245425            G2plt.PlotSASDPairDist(G2frame)
    5425             RefreshPlots(True)
    54265426            wx.CallAfter(UpdateModelsGrid,G2frame,data)
    54275427           
     
    54325432      A New Algroithm for the Reconstruction of Protein Molecular Envelopes
    54335433      from X-ray Solution Scattering Data,
    5434       J. Badger, Jour. of Appl. Chrystallogr. 2019, XX, xxx-xxx.
    5435       doi: xxx''',
     5434      J. Badger, Jour. of Appl. Chrystallogr. 2019, 52, xxx-xxx.
     5435      doi: 10.1107/S1600576719009774''',
    54365436      caption='Program Shapes',style=wx.ICON_INFORMATION)
    54375437            data['Pair']['Result'] = []       #clear old results (if any) for now
     
    56075607        def OnMooreTerms(event):
    56085608            data['Pair']['Moore'] = int(round(Limits[1][1]*data['Pair']['MaxRadius']/np.pi))-1
    5609             wx.CallAfter(UpdateModelsGrid,G2frame,data)                               
     5609            wx.CallAfter(UpdateModelsGrid,G2frame,data)
     5610
     5611        def OnNewVal(invalid,value,tc):
     5612            if invalid: return
     5613            parmDict = {'Rg':data['Pair']['MaxRadius']/2.5,'G':data['Pair']['Dist G'],
     5614                'B':data['Pair'].get('Dist B',Profile[1][-1]*Profile[0][-1]**4),
     5615                'Back':data['Back'][0]}
     5616            Profile[2] = G2sasd.getSASDRg(Profile[0],parmDict)
     5617            RefreshPlots(True)
    56105618           
    56115619        pairSizer = wx.BoxSizer(wx.VERTICAL)
     
    56205628        binSizer.Add(nbins,0,WACV)
    56215629        binSizer.Add(wx.StaticText(G2frame.dataWindow,label=' Max diam.: '),0,WACV)
    5622         maxdiam = G2G.ValidatedTxtCtrl(G2frame.dataWindow,data['Pair'],'MaxRadius',min=10.,nDig=(10,1))
     5630        maxdiam = G2G.ValidatedTxtCtrl(G2frame.dataWindow,data['Pair'],'MaxRadius',min=10.,nDig=(10,1),OnLeave=OnNewVal)
    56235631        binSizer.Add(maxdiam,0,WACV)
    56245632        maxest = wx.Button(G2frame.dataWindow,label='Make estimate')
     
    56285636        pairSizer.Add((5,5),0)
    56295637        fitSizer = wx.BoxSizer(wx.HORIZONTAL)
    5630         methods = ['Regularization','Moore']
     5638        methods = ['Moore',]            #'Regularization',
    56315639        fitSizer.Add(wx.StaticText(G2frame.dataWindow,label='Fitting method: '),0,WACV)
    56325640        method = wx.ComboBox(G2frame.dataWindow,value=data['Pair']['Method'],choices=methods,
     
    56345642        method.Bind(wx.EVT_COMBOBOX,OnMethod)
    56355643        fitSizer.Add(method,0,WACV)
     5644        if data['Pair']['Method'] == 'Moore':
     5645            fitSizer.Add(wx.StaticText(G2frame.dataWindow,label=" P.B. Moore, J. Appl. Cryst., 13, 168-175 (1980)"),0,WACV)
     5646        else:
     5647            fitSizer.Add(wx.StaticText(G2frame.dataWindow,label=" D.I. Svergun, J. Appl. Cryst., 24, 485-492 (1991)"),0,WACV)
    56365648        pairSizer.Add(fitSizer,0,WACV)
    56375649        if 'Moore' in data['Pair']['Method']:
    56385650            mooreSizer = wx.BoxSizer(wx.HORIZONTAL)
    56395651            mooreSizer.Add(wx.StaticText(G2frame.dataWindow,label='Number of functions: '),0,WACV)
    5640             moore = G2G.ValidatedTxtCtrl(G2frame.dataWindow,data['Pair'],'Moore',min=10)
     5652            moore = G2G.ValidatedTxtCtrl(G2frame.dataWindow,data['Pair'],'Moore',min=2)
    56415653            mooreSizer.Add(moore,0,WACV)
    56425654            mooreterms = wx.Button(G2frame.dataWindow,label = 'Auto determine?')
     
    56715683            selAtoms = Atoms[2*r+(c-1)]
    56725684            pattern = Patterns[r]
    5673             data['Pair']['Pair Calc'] = np.array([PRcalc[r][0],PRcalc[r][2]]).T
     5685            prCalc = PRcalc[r][2]
     5686            prDelt= np.diff(PRcalc[r][0])[0]
     5687            prsum = np.sum(prCalc)
     5688            prCalc /= prsum*prDelt
     5689            data['Pair']['Pair Calc'] = np.array([PRcalc[r][0],prCalc]).T
    56745690            print('%s %d'%('num. beads',len(selAtoms[1])))
    56755691            print('%s %.3f'%('selected r value',pattern[-1]))
    56765692            print('%s %.3f'%('selected Delta P(r)',PRcalc[r][-1]))
     5693            G2plt.PlotBeadModel(G2frame,selAtoms,plotDefaults)
     5694            RefreshPlots(True)
    56775695            G2plt.PlotSASDPairDist(G2frame)
    5678             RefreshPlots(True)
    56795696           
    5680             G2plt.PlotBeadModel(G2frame,selAtoms,plotDefaults)
    56815697           
    56825698           
     
    56965712        inflate = G2G.ValidatedTxtCtrl(G2frame.dataWindow,data['Shapes'],'inflateV',min=1.,max=1.4,nDig=(10,2))
    56975713        parmSizer.Add(inflate,0,WACV)
    5698         parmSizer.Add(wx.StaticText(G2frame.dataWindow,label=' Axial symmetry (1-6): '),0,WACV)       
    5699         symm = G2G.ValidatedTxtCtrl(G2frame.dataWindow,data['Shapes'],'Symm',min=1,max=6)
     5714        parmSizer.Add(wx.StaticText(G2frame.dataWindow,label=' Axial symmetry (1-12): '),0,WACV)       
     5715        symm = G2G.ValidatedTxtCtrl(G2frame.dataWindow,data['Shapes'],'Symm',min=1,max=12)
    57005716        parmSizer.Add(symm,0,WACV)
    57015717#3rd row
     
    60746090    backSizer.Add(wx.StaticText(G2frame.dataWindow,label=' Background:'),0,WACV)
    60756091    backVal = G2G.ValidatedTxtCtrl(G2frame.dataWindow,data['Back'],0,
    6076         nDig=(10,3),typeHint=float,OnLeave=OnBackChange)
     6092        nDig=(10,3,'g'),OnLeave=OnBackChange)
    60776093    backSizer.Add(backVal,0,WACV)
    60786094    backVar = wx.CheckBox(G2frame.dataWindow,label='Refine?')
  • trunk/GSASIIsasd.py

    r4077 r4094  
    666666    DISTANCE_LIMIT_FACTOR = 0.1                 # limitation on df to constrain runaways
    667667   
    668     MAX_MOVE_LOOPS    = 500                     # for no solution in routine: move,
     668    MAX_MOVE_LOOPS    = 5000                     # for no solution in routine: move,
    669669    MOVE_PASSES       = 0.001                   # convergence test in routine: move
    670670
     
    11281128        if ifBack:
    11291129            N += 1
    1130         MPV = np.zeros(N)
    1131         MPV[0] = Q[Ibeg]
     1130        MPV = np.ones(N)*1.e-5
     1131#        MPV[0] = Q[Ibeg]
    11321132        dmax = pairData['MaxRadius']
    11331133        result = so.leastsq(calcSASD,MPV,full_output=True,epsfcn=1.e-8,   #ftol=Ftol,
     
    11411141            Back = 0.
    11421142        chisq = np.sum(result[2]['fvec']**2)
    1143         Ic[Ibeg:Ifin] = MooreIOREFF(MPVR,Q[Ibeg:Ifin],dmax)+Ifb+Back
     1143        covM = result[1]
     1144        Ic[Ibeg:Ifin] = MooreIOREFF(MPVR,Q[Ibeg:Ifin],dmax)+Ifb[Ibeg:Ifin]+Back
    11441145        ncalc = result[2]['nfev']
    11451146        GOF = chisq/(Ifin-Ibeg-N)
     
    11481149        print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(ncalc,Ifin-Ibeg,N))
    11491150        print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF))
     1151        if len(covM):
     1152            sig = np.sqrt(np.diag(covM)*GOF)
     1153            for val,esd in zip(result[0],sig):
     1154                print(' parameter: %.4g esd: %.4g'%(val,esd))           
    11501155        BinMag = MoorePOR(MPVR,Bins,dmax)/2.
    11511156        return Bins,Dbins,BinMag
    11521157   
    1153     def CalcRegular():
    1154         Bins = np.linspace(0.,pairData['MaxRadius'],pairData['NBins']+1,True)
    1155         Dbins = np.diff(Bins)
    1156         Bins = Bins[:-1]+Dbins/2.
    1157 
    1158 #test       
    1159         midBin = pairData['MaxRadius']/4.
    1160         wid = midBin/4.
    1161         BinMag  = 1./np.sqrt(2.*np.pi*wid*2)*np.exp(-(midBin-Bins)**2/(2.*wid**2))
    1162         return Bins,Dbins,BinMag
    1163                
    11641158    pairData = data['Pair']
    11651159   
    1166     if pairData['Method'] == 'Regularization':
     1160    if pairData['Method'] == 'Regularization':      #not used
    11671161        print('Regularization calc; dummy Gaussian - TBD')
    1168         Bins,Dbins,BinMag = CalcRegular()
     1162        pairData['Method'] = 'Moore'
    11691163       
    11701164       
    11711165    elif pairData['Method'] == 'Moore':
    1172         Bins,Dbins,BinMag = CalcMoore()       
    1173    
    1174     data['Pair']['Distribution'] = [Bins,Dbins,BinMag/(2.*Dbins)]
     1166        Bins,Dbins,BinMag = CalcMoore()
     1167        BinSum = np.sum(BinMag*Dbins)
     1168        BinMag /= BinSum       
     1169   
     1170    data['Pair']['Distribution'] = [Bins,Dbins,BinMag]      #/(2.*Dbins)]
     1171    if 'Pair Calc' in data['Pair']:
     1172        del data['Pair']['Pair Calc']
    11751173   
    11761174   
     
    13901388        return False,0,0,0,0,0,0,Msg
    13911389
     1390def getSASDRg(Q,parmDict):
     1391    Ic = np.zeros_like(Q)
     1392    Rg,G,B = parmDict['Rg'],parmDict['G'],parmDict['B']
     1393    Qstar = Q/(scsp.erf(Q*Rg/np.sqrt(6)))**3
     1394    Guin = G*np.exp(-(Q*Rg)**2/3)
     1395    Porod = (B/Qstar**4)        #*np.exp(-(Q*B)**2/3)
     1396    Ic += Guin+Porod+parmDict['Back']
     1397    return Ic
     1398       
    13921399def RgFit(Profile,ProfDict,Limits,Sample,Model):
    13931400    print('unified fit single Rg to data to estimate a size')
     
    14031410        parmDict['Back'] = Back[0]
    14041411        pairData = Model['Pair']
    1405         parmDict['G'] = pairData.get('Dist G',100.)
     1412        parmDict['G'] = pairData.get('Dist G',Io[Ibeg])
    14061413        parmDict['Rg'] = pairData['MaxRadius']/2.5
    1407         varyList  += ['G','Rg',]
    1408         values += [parmDict['G'],parmDict['Rg'],]
     1414        parmDict['B'] = pairData.get('Dist B',Io[Ifin-1]*Q[Ifin-1]**4)
     1415        varyList  += ['G','Rg','B']
     1416        values += [parmDict['G'],parmDict['Rg'],parmDict['B']]
    14091417        return parmDict,varyList,values
    14101418
    14111419    def calcSASD(values,Q,Io,wt,Ifb,parmDict,varyList):
    14121420        parmDict.update(zip(varyList,values))
    1413         M = np.sqrt(wt)*(getSASD(Q,parmDict)+Ifb-Io)
     1421        M = np.sqrt(wt)*(getSASDRg(Q,parmDict)-Io)
    14141422        return M
    14151423   
    1416     def getSASD(Q,parmDict):
    1417         Ic = np.zeros_like(Q)
    1418         Rg,G = parmDict['Rg'],parmDict['G']
    1419         Guin = G*np.exp(-(Q*Rg)**2/3)
    1420         Ic += Guin
    1421         Ic += parmDict['Back']  #/parmDict['Scale']
    1422         return Ic
    1423        
    14241424    def SetModelParms():
    14251425        print (' Refined parameters: ')
     
    14301430        pairData['Dist G'] = parmDict['G']
    14311431        pairData['MaxRadius'] = parmDict['Rg']*2.5
     1432        pairData['Dist B'] = parmDict['B']
    14321433        for item in varyList:
    14331434            print (' %15s: %15.4g esd: %15.4g'%(item,parmDict[item],sigDict[item]))
     
    14411442    Ic[:] = 0
    14421443    parmDict,varyList,values = GetModelParms()
    1443     result = so.leastsq(calcSASD,values,full_output=True,epsfcn=1.e-8,   #ftol=Ftol,
     1444    result = so.leastsq(calcSASD,values,full_output=True,epsfcn=1.e-12,factor=0.1,  #ftol=Ftol,
    14441445        args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Ifb[Ibeg:Ifin],parmDict,varyList))
    1445     parmDict.update(zip(varyList,result[0]))
     1446    parmDict.update(dict(zip(varyList,result[0])))
    14461447    chisq = np.sum(result[2]['fvec']**2)
    14471448    ncalc = result[2]['nfev']
     
    14501451    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
    14511452    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
    1452     Ic[Ibeg:Ifin] = getSASD(Q[Ibeg:Ifin],parmDict)
     1453    Ic[Ibeg:Ifin] = getSASDRg(Q[Ibeg:Ifin],parmDict)
    14531454    Msg = 'Failed to converge'
    14541455    try:
  • trunk/Substances.py

    r2776 r4094  
    5757'SrTiO3':{'Elements':{'Sr':{'Num':1},'Ti':{'Num':1},'O':{'Num':1}},'Volume':26.71},
    5858'V':{'Elements':{'V':{'Num':1}},'Volume':19.26},
     59'protein':{'Elements':{'C':{'Num':9.25},'N':{'Num':2.34},'O':{'Num':2.77},'S':{'Num':0.22},'H':{'Num':14.3}},'Volume':288.8}, #insulin - typical?
    5960}
    6061# they should not be duplicated in the UserSubstances.py file:
Note: See TracChangeset for help on using the changeset viewer.