source: trunk/GSASIIsasd.py @ 1269

Last change on this file since 1269 was 1269, checked in by vondreele, 8 years ago

more model fitting in SASD
add a possible for plotting hkls from powder data - commented out fro now

File size: 40.5 KB
Line 
1#/usr/bin/env python
2# -*- coding: utf-8 -*-
3'''
4*GSASII small angle calculation module*
5==================================
6
7'''
8########### SVN repository information ###################
9# $Date: 2014-01-09 11:09:53 -0600 (Thu, 09 Jan 2014) $
10# $Author: vondreele $
11# $Revision: 1186 $
12# $URL: https://subversion.xray.aps.anl.gov/pyGSAS/trunk/GSASIIsasd.py $
13# $Id: GSASIIsasd.py 1186 2014-01-09 17:09:53Z vondreele $
14########### SVN repository information ###################
15import os
16import sys
17import math
18import time
19
20import numpy as np
21import scipy as sp
22import numpy.linalg as nl
23from numpy.fft import ifft, fft, fftshift
24import scipy.special as scsp
25import scipy.interpolate as si
26import scipy.stats as st
27import scipy.optimize as so
28
29import GSASIIpath
30GSASIIpath.SetVersionNumber("$Revision: 1186 $")
31import GSASIIlattice as G2lat
32import GSASIIspc as G2spc
33import GSASIIElem as G2elem
34import GSASIIgrid as G2gd
35import GSASIIIO as G2IO
36import GSASIImath as G2mth
37import pypowder as pyd
38
39# trig functions in degrees
40sind = lambda x: math.sin(x*math.pi/180.)
41asind = lambda x: 180.*math.asin(x)/math.pi
42tand = lambda x: math.tan(x*math.pi/180.)
43atand = lambda x: 180.*math.atan(x)/math.pi
44atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
45cosd = lambda x: math.cos(x*math.pi/180.)
46acosd = lambda x: 180.*math.acos(x)/math.pi
47rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
48#numpy versions
49npsind = lambda x: np.sin(x*np.pi/180.)
50npasind = lambda x: 180.*np.arcsin(x)/math.pi
51npcosd = lambda x: np.cos(x*math.pi/180.)
52npacosd = lambda x: 180.*np.arccos(x)/math.pi
53nptand = lambda x: np.tan(x*math.pi/180.)
54npatand = lambda x: 180.*np.arctan(x)/np.pi
55npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
56npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave
57npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
58   
59###############################################################################
60#### Particle form factors
61###############################################################################
62
63def SphereFF(Q,R,args=()):
64    ''' 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)
67    param array args: ignored
68    returns float: form factors as array as needed
69    '''
70    QR = Q[:,np.newaxis]*R
71    return (3./(QR**3))*(np.sin(QR)-(QR*np.cos(QR)))
72   
73def SpheroidFF(Q,R,args):
74    ''' Compute form factor of cylindrically symmetric ellipsoid (spheroid)
75    - can use numpy arrays for R & AR; will return corresponding numpy array
76    param float Q : Q value array (usually in A-1)
77    param float R: radius along 2 axes of spheroid
78    param array args: [float AR]: aspect ratio so 3rd axis = R*AR
79    returns float: form factors as array as needed
80    '''
81    NP = 50
82    AR = args[0]
83    if 0.99 < AR < 1.01:
84        return SphereFF(Q,R,0)
85    else:
86        cth = np.linspace(0,1.,NP)
87        Rct = R[:,np.newaxis]*np.sqrt(1.+(AR**2-1.)*cth**2)
88        return np.sqrt(np.sum(SphereFF(Q[:,np.newaxis],Rct,0)**2,axis=2)/NP)
89           
90def CylinderFF(Q,R,args):
91    ''' Compute form factor for cylinders - can use numpy arrays
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)
95    returns float: form factor
96    '''
97    L = args[0]
98    NP = 200
99    alp = np.linspace(0,np.pi/2.,NP)
100    QL = Q[:,np.newaxis]*L
101    LBessArg = 0.5*QL[:,:,np.newaxis]**np.cos(alp)
102    LBess = np.where(LBessArg<1.e-6,1.,np.sin(LBessArg)/LBessArg)
103    QR = Q[:,np.newaxis]*R
104    SBessArg = QR[:,:,np.newaxis]*np.sin(alp)
105    SBess = np.where(SBessArg<1.e-6,0.5,scsp.jv(1,SBessArg)/SBessArg)
106    LSBess = LBess*SBess
107    return np.sqrt(2.*np.pi*np.sum(np.sin(alp)*LSBess**2,axis=2)/NP)
108   
109def CylinderDFF(Q,L,args):
110    ''' Compute form factor for cylinders - can use numpy arrays
111    param float Q: Q value array (A-1)
112    param float L: cylinder half length (A)
113    param array args: [float R]: cylinder radius (A)
114    returns float: form factor
115    '''
116    R = args[0]
117    return CylinderFF(Q,R,args=[2.*L])   
118   
119def CylinderARFF(Q,R,args): 
120    ''' Compute form factor for cylinders - can use numpy arrays
121    param float Q: Q value array (A-1)
122    param float R: cylinder radius (A)
123    param array args: [float AR]: cylinder aspect ratio = L/D = L/2R
124    returns float: form factor
125    '''
126    AR = args[0]
127    return CylinderFF(Q,R,args=[2.*R*AR])   
128   
129def UniSphereFF(Q,R,args=0):
130    ''' Compute form factor for unified sphere - can use numpy arrays
131    param float Q: Q value array (A-1)
132    param float R: cylinder radius (A)
133    param array args: ignored
134    returns float: form factor
135    '''
136    Rg = np.sqrt(3./5.)*R
137    B = np.pi*1.62/(Rg**4)    #are we missing *np.pi? 1.62 = 6*(3/5)**2/(4/3) sense?
138    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg/np.sqrt(6)))**3
139    return np.sqrt(np.exp((-Q[:,np.newaxis]**2*Rg**2)/3.)+(B/QstV**4))
140   
141def UniRodFF(Q,R,args):
142    ''' Compute form factor for unified rod - can use numpy arrays
143    param float Q: Q value array (A-1)
144    param float R: cylinder radius (A)
145    param array args: [float R]: cylinder radius (A)
146    returns float: form factor
147    '''
148    L = args[0]
149    Rg2 = np.sqrt(R**2/2+L**2/12)
150    B2 = np.pi/L
151    Rg1 = np.sqrt(3.)*R/2.
152    G1 = (2./3.)*R/L
153    B1 = 4.*(L+R)/(R**3*L**2)
154    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg2/np.sqrt(6)))**3
155    FF = np.exp(-Q[:,np.newaxis]**2*Rg2**2/3.)+(B2/QstV)*np.exp(-Rg1**2*Q[:,np.newaxis]**2/3.)
156    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg1/np.sqrt(6)))**3
157    FF += G1*np.exp(-Q[:,np.newaxis]**2*Rg1**2/3.)+(B1/QstV**4)
158    return np.sqrt(FF)
159   
160def UniRodARFF(Q,R,args):
161    ''' Compute form factor for unified rod of fixed aspect ratio - can use numpy arrays
162    param float Q: Q value array (A-1)
163    param float R: cylinder radius (A)
164    param array args: [float AR]: cylinder aspect ratio = L/D = L/2R
165    returns float: form factor
166    '''
167    AR = args[0]
168    return UniRodFF(Q,R,args=[AR*R,])
169   
170def UniDiskFF(Q,R,args):
171    ''' Compute form factor for unified disk - can use numpy arrays
172    param float Q: Q value array (A-1)
173    param float R: cylinder radius (A)
174    param array args: [float T]: disk thickness (A)
175    returns float: form factor
176    '''
177    T = args[0]
178    Rg2 = np.sqrt(R**2/2.+T**2/12.)
179    B2 = 2./R**2
180    Rg1 = np.sqrt(3.)*T/2.
181    RgC2 = 1.1*Rg1
182    G1 = (2./3.)*(T/R)**2
183    B1 = 4.*(T+R)/(R**3*T**2)
184    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg2/np.sqrt(6)))**3
185    FF = np.exp(-Q[:,np.newaxis]**2*Rg2**2/3.)+(B2/QstV**2)*np.exp(-RgC2**2*Q[:,np.newaxis]**2/3.)
186    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg1/np.sqrt(6)))**3
187    FF += G1*np.exp(-Q[:,np.newaxis]**2*Rg1**2/3.)+(B1/QstV**4)
188    return np.sqrt(FF)
189   
190def UniTubeFF(Q,R,args):
191    ''' Compute form factor for unified disk - can use numpy arrays
192    param float Q: Q value array (A-1)
193    param float R: cylinder radius (A)
194    param array args: [float L,T]: tube length & wall thickness(A)
195    returns float: form factor
196    '''
197    L,T = args[:2]
198    Ri = R-T
199    DR2 = R**2-Ri**2
200    Vt = np.pi*DR2*L
201    Rg3 = np.sqrt(DR2/2.+L**2/12.)
202    B1 = 4.*np.pi**2*(DR2+L*(R+Ri))/Vt**2
203    B2 = np.pi**2*T/Vt
204    B3 = np.pi/L
205    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg3/np.sqrt(6)))**3
206    FF = np.exp(-Q[:,np.newaxis]**2*Rg3**2/3.)+(B3/QstV)*np.exp(-Q[:,np.newaxis]**2*R**2/3.)
207    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*R/np.sqrt(6)))**3
208    FF += (B2/QstV**2)*np.exp(-Q[:,np.newaxis]**2*T**2/3.)
209    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*T/np.sqrt(6)))**3
210    FF += B1/QstV**4
211    return np.sqrt(FF)
212
213###############################################################################
214#### Particle volumes
215###############################################################################
216
217def SphereVol(R,args=()):
218    ''' Compute volume of sphere
219    - numpy array friendly
220    param float R: sphere radius
221    param array args: ignored
222    returns float: volume
223    '''
224    return (4./3.)*np.pi*R**3
225
226def SpheroidVol(R,args):
227    ''' Compute volume of cylindrically symmetric ellipsoid (spheroid)
228    - numpy array friendly
229    param float R: radius along 2 axes of spheroid
230    param array args: [float AR]: aspect ratio so radius of 3rd axis = R*AR
231    returns float: volume
232    '''
233    AR = args[0]
234    return AR*SphereVol(R)
235   
236def CylinderVol(R,args):
237    ''' Compute cylinder volume for radius & length
238    - numpy array friendly
239    param float R: diameter (A)
240    param array args: [float L]: length (A)
241    returns float:volume (A^3)
242    '''
243    L = args[0]
244    return np.pi*L*R**2
245   
246def CylinderDVol(L,args):
247    ''' Compute cylinder volume for length & diameter
248    - numpy array friendly
249    param float: L half length (A)
250    param array args: [float D]: diameter (A)
251    returns float:volume (A^3)
252    '''
253    D = args[0]
254    return CylinderVol(D/2.,[2.*L,])
255   
256def CylinderARVol(R,args):
257    ''' Compute cylinder volume for radius & aspect ratio = L/D
258    - numpy array friendly
259    param float: R radius (A)
260    param array args: [float AR]: =L/D=L/2R aspect ratio
261    returns float:volume
262    '''
263    AR = args[0]
264    return CylinderVol(R,[2.*R*AR,])
265   
266def UniSphereVol(R,args=()):
267    ''' Compute volume of sphere
268    - numpy array friendly
269    param float R: sphere radius
270    param array args: ignored
271    returns float: volume
272    '''
273    return SphereVol(R)
274   
275def UniRodVol(R,args):
276    ''' Compute cylinder volume for radius & length
277    - numpy array friendly
278    param float R: diameter (A)
279    param array args: [float L]: length (A)
280    returns float:volume (A^3)
281    '''
282    L = args[0]
283    return CylinderVol(R,[L,])
284   
285def UniRodARVol(R,args):
286    ''' Compute rod volume for radius & aspect ratio
287    - numpy array friendly
288    param float R: diameter (A)
289    param array args: [float AR]: =L/D=L/2R aspect ratio
290    returns float:volume (A^3)
291    '''
292    AR = args[0]
293    return CylinderARVol(R,[AR,])
294   
295def UniDiskVol(R,args):
296    ''' Compute disk volume for radius & thickness
297    - numpy array friendly
298    param float R: diameter (A)
299    param array args: [float T]: thickness
300    returns float:volume (A^3)
301    '''
302    T = args[0]
303    return CylinderVol(R,[T,])
304   
305def UniTubeVol(R,args):
306    ''' Compute tube volume for radius, length & wall thickness
307    - numpy array friendly
308    param float R: diameter (A)
309    param array args: [float L,T]: tube length & wall thickness(A)
310    returns float: volume (A^3) of tube wall
311    '''
312    L,T = arg[:2]
313    return CylinderVol(R,[L,])-CylinderVol(R-T,[L,])
314   
315################################################################################
316#### Distribution functions & their cumulative fxns
317################################################################################
318
319def LogNormalDist(x,pos,args):
320    ''' Standard LogNormal distribution - numpy friendly on x axis
321    ref: http://www.itl.nist.gov/div898/handbook/index.htm 1.3.6.6.9
322    param float x: independent axis (can be numpy array)
323    param float pos: location of distribution
324    param float scale: width of distribution (m)
325    param float shape: shape - (sigma of log(LogNormal))
326    returns float: LogNormal distribution
327    '''
328    scale,shape = args
329    return np.exp(-np.log((x-pos)/scale)**2/(2.*shape**2))/(np.sqrt(2.*np.pi)*(x-pos)*shape)
330   
331def GaussDist(x,pos,args):
332    ''' Standard Normal distribution - numpy friendly on x axis
333    param float x: independent axis (can be numpy array)
334    param float pos: location of distribution
335    param float scale: width of distribution (sigma)
336    param float shape: not used
337    returns float: Normal distribution
338    '''
339    scale = args[0]
340    return (1./(scale*np.sqrt(2.*np.pi)))*np.exp(-(x-pos)**2/(2.*scale**2))
341   
342def LSWDist(x,pos,args=[]):
343    ''' Lifshitz-Slyozov-Wagner Ostwald ripening distribution - numpy friendly on x axis
344    ref:
345    param float x: independent axis (can be numpy array)
346    param float pos: location of distribution
347    param float scale: not used
348    param float shape: not used
349    returns float: LSW distribution   
350    '''
351    redX = x/pos
352    result = (81.*2**(-5/3.))*(redX**2*np.exp(-redX/(1.5-redX)))/((1.5-redX)**(11/3.)*(3.-redX)**(7/3.))
353   
354    return np.nan_to_num(result/pos)
355   
356def SchulzZimmDist(x,pos,args):
357    ''' Schulz-Zimm macromolecule distribution - numpy friendly on x axis
358    ref: http://goldbook.iupac.org/S05502.html
359    param float x: independent axis (can be numpy array)
360    param float pos: location of distribution
361    param float scale: width of distribution (sigma)
362    param float shape: not used
363    returns float: Schulz-Zimm distribution
364    '''
365    scale = args[0]
366    b = (2.*pos/scale)**2
367    a = b/pos
368    if b<70:    #why bother?
369        return (a**(b+1.))*x**b*np.exp(-a*x)/scsp.gamma(b+1.)
370    else:
371        return np.exp((b+1.)*np.log(a)-scsp.gammaln(b+1.)+b*np.log(x)-(a*x))
372           
373def LogNormalCume(x,pos,args):
374    ''' Standard LogNormal cumulative distribution - numpy friendly on x axis
375    ref: http://www.itl.nist.gov/div898/handbook/index.htm 1.3.6.6.9
376    param float x: independent axis (can be numpy array)
377    param float pos: location of distribution
378    param float scale: width of distribution (sigma)
379    param float shape: shape parameter
380    returns float: LogNormal cumulative distribution
381    '''
382    scale,shape = args
383    lredX = np.log((x-pos)/scale)
384    return (scsp.erf((lredX/shape)/np.sqrt(2.))+1.)/2.
385   
386def GaussCume(x,pos,args):
387    ''' Standard Normal cumulative distribution - numpy friendly on x axis
388    param float x: independent axis (can be numpy array)
389    param float pos: location of distribution
390    param float scale: width of distribution (sigma)
391    param float shape: not used
392    returns float: Normal cumulative distribution
393    '''
394    scale = args[0]
395    redX = (x-pos)/scale
396    return (scsp.erf(redX/np.sqrt(2.))+1.)/2.
397   
398def LSWCume(x,pos,args=[]):
399    ''' Lifshitz-Slyozov-Wagner Ostwald ripening cumulative distribution - numpy friendly on x axis
400    param float x: independent axis (can be numpy array)
401    param float pos: location of distribution
402    param float scale: not used
403    param float shape: not used
404    returns float: LSW cumulative distribution
405    '''
406    nP = 500
407    xMin,xMax = [0.,2.*pos]
408    X = np.linspace(xMin,xMax,nP,True)
409    fxn = LSWDist(X,pos)
410    mat = np.outer(np.ones(nP),fxn)
411    cume = np.sum(np.tril(mat),axis=1)/np.sum(fxn)
412    return np.interp(x,X,cume,0,1)
413   
414def SchulzZimmCume(x,pos,args):
415    ''' Schulz-Zimm cumulative distribution - numpy friendly on x axis
416    param float x: independent axis (can be numpy array)
417    param float pos: location of distribution
418    param float scale: width of distribution (sigma)
419    param float shape: not used
420    returns float: Normal distribution
421    '''
422    scale = args[0]
423    nP = 500
424    xMin = np.max([0.,pos-4.*scale])
425    xMax = np.min([pos+4.*scale,1.e5])
426    X = np.linspace(xMin,xMax,nP,True)
427    fxn = LSWDist(X,pos)
428    mat = np.outer(np.ones(nP),fxn)
429    cume = np.sum(np.tril(mat),axis=1)/np.sum(fxn)
430    return np.interp(x,X,cume,0,1)
431   
432    return []
433   
434         
435################################################################################
436##### SB-MaxEnt
437################################################################################
438
439def G_matrix(q,r,contrast,FFfxn,Volfxn,args=()):
440    '''Calculates the response matrix :math:`G(Q,r)`
441   
442    :param float q: :math:`Q`
443    :param float r: :math:`r`
444    :param float contrast: :math:`|\\Delta\\rho|^2`, the scattering contrast
445    :param function FFfxn: form factor function FF(q,r,args)
446    :param function Volfxn: volume function Vol(r,args)
447    :returns float: G(Q,r)
448    '''
449    FF = FFfxn(q,r,args)
450    Vol = Volfxn(r,args)
451    return 1.e-4*(contrast*Vol*FF**2).T     #10^-20 vs 10^-24
452   
453'''
454sbmaxent
455
456Entropy maximization routine as described in the article
457J Skilling and RK Bryan; MNRAS 211 (1984) 111 - 124.
458("MNRAS": "Monthly Notices of the Royal Astronomical Society")
459
460:license: Copyright (c) 2013, UChicago Argonne, LLC
461:license: This file is distributed subject to a Software License Agreement found
462     in the file LICENSE that is included with this distribution.
463
464References:
465
4661. J Skilling and RK Bryan; MON NOT R ASTR SOC 211 (1984) 111 - 124.
4672. JA Potton, GJ Daniell, and BD Rainford; Proc. Workshop
468   Neutron Scattering Data Analysis, Rutherford
469   Appleton Laboratory, UK, 1986; ed. MW Johnson,
470   IOP Conference Series 81 (1986) 81 - 86, Institute
471   of Physics, Bristol, UK.
4723. ID Culverwell and GP Clarke; Ibid. 87 - 96.
4734. JA Potton, GK Daniell, & BD Rainford,
474   J APPL CRYST 21 (1988) 663 - 668.
4755. JA Potton, GJ Daniell, & BD Rainford,
476   J APPL CRYST 21 (1988) 891 - 897.
477
478'''
479
480class MaxEntException(Exception): 
481    '''Any exception from this module'''
482    pass
483
484def MaxEnt_SB(datum, sigma, G, base, IterMax, image_to_data=None, data_to_image=None, report=False):
485    '''
486    do the complete Maximum Entropy algorithm of Skilling and Bryan
487   
488    :param float datum[]:
489    :param float sigma[]:
490    :param float[][] G: transformation matrix
491    :param float base[]:
492    :param int IterMax:
493    :param obj image_to_data: opus function (defaults to opus)
494    :param obj data_to_image: tropus function (defaults to tropus)
495   
496    :returns float[]: :math:`f(r) dr`
497    '''
498       
499    TEST_LIMIT        = 0.05                    # for convergence
500    CHI_SQR_LIMIT     = 0.01                    # maximum difference in ChiSqr for a solution
501    SEARCH_DIRECTIONS = 3                       # <10.  This code requires value = 3
502    RESET_STRAYS      = 1                       # was 0.001, correction of stray negative values
503    DISTANCE_LIMIT_FACTOR = 0.1                 # limitation on df to constrain runaways
504   
505    MAX_MOVE_LOOPS    = 500                     # for no solution in routine: move,
506    MOVE_PASSES       = 0.001                   # convergence test in routine: move
507
508    def tropus (data, G):
509        '''
510        tropus: transform data-space -> solution-space:  [G] * data
511       
512        default definition, caller can use this definition or provide an alternative
513       
514        :param float[M] data: observations, ndarray of shape (M)
515        :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
516        :returns float[N]: calculated image, ndarray of shape (N)
517        '''
518        return G.dot(data)
519
520    def opus (image, G):
521        '''
522        opus: transform solution-space -> data-space:  [G]^tr * image
523       
524        default definition, caller can use this definition or provide an alternative
525       
526        :param float[N] image: solution, ndarray of shape (N)
527        :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
528        :returns float[M]: calculated data, ndarray of shape (M)
529        '''
530        return np.dot(G.T,image)    #G.transpose().dot(image)
531
532    def Dist(s2, beta):
533        '''measure the distance of this possible solution'''
534        w = 0
535        n = beta.shape[0]
536        for k in range(n):
537            z = -sum(s2[k] * beta)
538            w += beta[k] * z
539        return w
540   
541    def ChiNow(ax, c1, c2, s1, s2):
542        '''
543        ChiNow
544       
545        :returns tuple: (ChiNow computation of ``w``, beta)
546        '''
547       
548        bx = 1 - ax
549        a =   bx * c2  -  ax * s2
550        b = -(bx * c1  -  ax * s1)
551   
552        beta = ChoSol(a, b)
553        w = 1.0
554        for k in range(SEARCH_DIRECTIONS):
555            w += beta[k] * (c1[k] + 0.5*sum(c2[k] * beta))
556        return w, beta
557   
558    def ChoSol(a, b):
559        '''
560        ChoSol: ? chop the solution vectors ?
561       
562        :returns: new vector beta
563        '''
564        n = b.shape[0]
565        fl = np.zeros((n,n))
566        bl = np.zeros_like(b)
567       
568        #print_arr("ChoSol: a", a)
569        #print_vec("ChoSol: b", b)
570   
571        if (a[0][0] <= 0):
572            msg = "ChoSol: a[0][0] = " 
573            msg += str(a[0][0])
574            msg += '  Value must be positive'
575            raise MaxEntException(msg)
576   
577        # first, compute fl from a
578        # note fl is a lower triangular matrix
579        fl[0][0] = math.sqrt (a[0][0])
580        for i in (1, 2):
581            fl[i][0] = a[i][0] / fl[0][0]
582            for j in range(1, i+1):
583                z = 0.0
584                for k in range(j):
585                    z += fl[i][k] * fl[j][k]
586                    #print "ChoSol: %d %d %d  z = %lg" % ( i, j, k, z)
587                z = a[i][j] - z
588                if j == i:
589                    y = math.sqrt(z)
590                else:
591                    y = z / fl[j][j]
592                fl[i][j] = y
593        #print_arr("ChoSol: fl", fl)
594   
595        # next, compute bl from fl and b
596        bl[0] = b[0] / fl[0][0]
597        for i in (1, 2):
598            z = 0.0
599            for k in range(i):
600                z += fl[i][k] * bl[k]
601                #print "\t", i, k, z
602            bl[i] = (b[i] - z) / fl[i][i]
603        #print_vec("ChoSol: bl", bl)
604   
605        # last, compute beta from bl and fl
606        beta = np.empty((n))
607        beta[-1] = bl[-1] / fl[-1][-1]
608        for i in (1, 0):
609            z = 0.0
610            for k in range(i+1, n):
611                z += fl[k][i] * beta[k]
612                #print "\t\t", i, k, 'z=', z
613            beta[i] = (bl[i] - z) / fl[i][i]
614        #print_vec("ChoSol: beta", beta)
615   
616        return beta
617
618    def MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2):
619        '''
620        move beta one step closer towards the solution
621        '''
622        a_lower, a_upper = 0., 1.          # bracket  "a"
623        cmin, beta = ChiNow (a_lower, c1, c2, s1, s2)
624        #print "MaxEntMove: cmin = %g" % cmin
625        if cmin*chisq > chizer:
626            ctarg = (1.0 + cmin)/2
627        else:
628            ctarg = chizer/chisq
629        f_lower = cmin - ctarg
630        c_upper, beta = ChiNow (a_upper, c1, c2, s1, s2)
631        f_upper = c_upper - ctarg
632   
633        fx = 2*MOVE_PASSES      # just to start off
634        loop = 1
635        while abs(fx) >= MOVE_PASSES and loop <= MAX_MOVE_LOOPS:
636            a_new = (a_lower + a_upper) * 0.5           # search by bisection
637            c_new, beta = ChiNow (a_new, c1, c2, s1, s2)
638            fx = c_new - ctarg
639            # tighten the search range for the next pass
640            if f_lower*fx > 0:
641                a_lower, f_lower = a_new, fx
642            if f_upper*fx > 0:
643                a_upper, f_upper = a_new, fx
644            loop += 1
645   
646        if abs(fx) >= MOVE_PASSES or loop > MAX_MOVE_LOOPS:
647            msg = "MaxEntMove: Loop counter = " 
648            msg += str(MAX_MOVE_LOOPS)
649            msg += '  No convergence in alpha chop'
650            raise MaxEntException(msg)
651   
652        w = Dist (s2, beta);
653        m = SEARCH_DIRECTIONS
654        if (w > DISTANCE_LIMIT_FACTOR*fSum/blank):        # invoke the distance penalty, SB eq. 17
655            for k in range(m):
656                beta[k] *= math.sqrt (fSum/(blank*w))
657        chtarg = ctarg * chisq
658        return w, chtarg, loop, a_new, fx, beta
659       
660#MaxEnt_SB starts here   
661
662    if image_to_data == None:
663        image_to_data = opus
664    if data_to_image == None:
665        data_to_image = tropus
666    n   = len(base)
667    npt = len(datum)
668
669    # Note that the order of subscripts for
670    # "xi" and "eta" has been reversed from
671    # the convention used in the FORTRAN version
672    # to enable parts of them to be passed as
673    # as vectors to "image_to_data" and "data_to_image".
674    xi      = np.zeros((SEARCH_DIRECTIONS, n))
675    eta     = np.zeros((SEARCH_DIRECTIONS, npt))
676    beta    = np.zeros((SEARCH_DIRECTIONS))
677    # s1      = np.zeros((SEARCH_DIRECTIONS))
678    # c1      = np.zeros((SEARCH_DIRECTIONS))
679    s2      = np.zeros((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
680    c2      = np.zeros((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
681
682    # TODO: replace blank (scalar) with base (vector)
683    blank = sum(base) / len(base)   # use the average value of base
684
685    chizer, chtarg = npt*1.0, npt*1.0
686    f = base * 1.0                              # starting distribution is base
687
688    fSum  = sum(f)                              # find the sum of the f-vector
689    z = (datum - image_to_data (f, G)) / sigma  # standardized residuals, SB eq. 3
690    chisq = sum(z*z)                            # Chi^2, SB eq. 4
691
692    for iter in range(IterMax):
693        ox = -2 * z / sigma                        # gradient of Chi^2
694
695        cgrad = data_to_image (ox, G)              # cgrad[i] = del(C)/del(f[i]), SB eq. 8
696        sgrad = -np.log(f/base) / (blank*math.exp (1.0))  # sgrad[i] = del(S)/del(f[i])
697        snorm = math.sqrt(sum(f * sgrad*sgrad))    # entropy term, SB eq. 22
698        cnorm = math.sqrt(sum(f * cgrad*cgrad))    # ChiSqr term, SB eq. 22
699        tnorm = sum(f * sgrad * cgrad)             # norm for gradient term TEST
700
701        a = 1.0
702        b = 1.0 / cnorm
703        if iter == 0:
704            test = 0.0     # mismatch between entropy and ChiSquared gradients
705        else:
706            test = math.sqrt( ( 1.0 - tnorm/(snorm*cnorm) )/2 ) # SB eq. 37?
707            a = 0.5 / (snorm * test)
708            b *= 0.5 / test
709        xi[0] = f * cgrad / cnorm
710        xi[1] = f * (a * sgrad - b * cgrad)
711
712        eta[0] = image_to_data (xi[0], G);          # image --> data
713        eta[1] = image_to_data (xi[1], G);          # image --> data
714        ox = eta[1] / (sigma * sigma)
715        xi[2] = data_to_image (ox, G);              # data --> image
716        a = 1.0 / math.sqrt(sum(f * xi[2]*xi[2]))
717        xi[2] = f * xi[2] * a
718        eta[2] = image_to_data (xi[2], G)           # image --> data
719       
720#         print_arr("MaxEnt: eta.transpose()", eta.transpose())
721#         print_arr("MaxEnt: xi.transpose()", xi.transpose())
722
723        # prepare the search directions for the conjugate gradient technique
724        c1 = xi.dot(cgrad) / chisq                          # C_mu, SB eq. 24
725        s1 = xi.dot(sgrad)                          # S_mu, SB eq. 24
726#         print_vec("MaxEnt: c1", c1)
727#         print_vec("MaxEnt: s1", s1)
728
729        for k in range(SEARCH_DIRECTIONS):
730            for l in range(k+1):
731                c2[k][l] = 2 * sum(eta[k] * eta[l] / sigma/sigma) / chisq
732                s2[k][l] = -sum(xi[k] * xi[l] / f) / blank
733#         print_arr("MaxEnt: c2", c2)
734#         print_arr("MaxEnt: s2", s2)
735
736        # reflect across the body diagonal
737        for k, l in ((0,1), (0,2), (1,2)):
738            c2[k][l] = c2[l][k]                     #  M_(mu,nu)
739            s2[k][l] = s2[l][k]                     #  g_(mu,nu)
740 
741        beta[0] = -0.5 * c1[0] / c2[0][0]
742        beta[1] = 0.0
743        beta[2] = 0.0
744        if (iter > 0):
745            w, chtarg, loop, a_new, fx, beta = MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2)
746
747        f_old = f.copy()    # preserve the last image
748        f += xi.transpose().dot(beta)   # move the image towards the solution, SB eq. 25
749       
750        # As mentioned at the top of p.119,
751        # need to protect against stray negative values.
752        # In this case, set them to RESET_STRAYS * base[i]
753        #f = f.clip(RESET_STRAYS * blank, f.max())
754        for i in range(n):
755            if f[i] <= 0.0:
756                f[i] = RESET_STRAYS * base[i]
757        df = f - f_old
758        fSum = sum(f)
759        fChange = sum(df)
760
761        # calculate the normalized entropy
762        S = sum((f/fSum) * np.log(f/fSum))      # normalized entropy, S&B eq. 1
763        z = (datum - image_to_data (f, G)) / sigma  # standardized residuals
764        chisq = sum(z*z)                            # report this ChiSq
765
766        if report:
767            print " MaxEnt trial/max: %3d/%3d" % ((iter+1), IterMax)
768            print " Residual: %5.2lf%% Entropy: %8lg" % (100*test, S)
769            if iter > 0:
770                value = 100*( math.sqrt(chisq/chtarg)-1)
771            else:
772                value = 0
773      #      print " %12.5lg %10.4lf" % ( math.sqrt(chtarg/npt), value )
774            print " Function sum: %.6lg Change from last: %.2lf%%\n" % (fSum,100*fChange/fSum)
775
776        # See if we have finished our task.
777        # do the hardest test first
778        if (abs(chisq/chizer-1.0) < CHI_SQR_LIMIT) and  (test < TEST_LIMIT):
779            print ' Convergence achieved.'
780            return chisq,f,image_to_data(f, G)     # solution FOUND returns here
781    print ' No convergence! Try increasing Error multiplier.'
782    return chisq,f,image_to_data(f, G)       # no solution after IterMax iterations
783
784   
785###############################################################################
786#### IPG/TNNLS Routines
787###############################################################################
788
789def IPG(datum,sigma,G,Bins,Dbins,IterMax,Qvec=[],approach=0.8,Power=-1,report=False):
790    ''' An implementation of the Interior-Point Gradient method of
791    Michael Merritt & Yin Zhang, Technical Report TR04-08, Dept. of Comp. and
792    Appl. Math., Rice Univ., Houston, Texas 77005, U.S.A. found on the web at
793    http://www.caam.rice.edu/caam/trs/2004/TR04-08.pdf
794    Problem addressed: Total Non-Negative Least Squares (TNNLS)
795    :param float datum[]:
796    :param float sigma[]:
797    :param float[][] G: transformation matrix
798    :param int IterMax:
799    :param float Qvec: data positions for Power = 0-4
800    :param float approach: 0.8 default fitting parameter
801    :param int Power: 0-4 for Q^Power weighting, -1 to use input sigma
802   
803    '''
804    if Power < 0: 
805        GmatE = G/sigma[:np.newaxis]
806        IntE = datum/sigma
807        pwr = 0
808        QvecP = np.ones_like(datum)
809    else:
810        GmatE = G[:]
811        IntE = datum[:]
812        pwr = Power
813        QvecP = Qvec**pwr
814    Amat = GmatE*QvecP[:np.newaxis]
815    AAmat = np.inner(Amat,Amat)
816    Bvec = IntE*QvecP
817    Xw = np.ones_like(Bins)*1.e-32
818    calc = np.dot(G.T,Xw)
819    nIter = 0
820    err = 10.
821    while (nIter<IterMax) and (err > 1.):
822        #Step 1 in M&Z paper:
823        Qk = np.inner(AAmat,Xw)-np.inner(Amat,Bvec)
824        Dk = Xw/np.inner(AAmat,Xw)
825        Pk = -Dk*Qk
826        #Step 2 in M&Z paper:
827        alpSt = -np.inner(Pk,Qk)/np.inner(Pk,np.inner(AAmat,Pk))
828        alpWv = -Xw/Pk
829        AkSt = [np.where(Pk[i]<0,np.min((approach*alpWv[i],alpSt)),alpSt) for i in range(Pk.shape[0])]
830        #Step 3 in M&Z paper:
831        shift = AkSt*Pk
832        Xw += shift
833        #done IPG; now results
834        nIter += 1
835        calc = np.dot(G.T,Xw)
836        chisq = np.sum(((datum-calc)/sigma)**2)
837        err = chisq/len(datum)
838        if report:
839            print ' Iteration: %d, chisq: %.3g, sum(shift^2): %.3g'%(nIter,chisq,np.sum(shift**2))
840    return chisq,Xw,calc
841
842###############################################################################
843#### SASD Utilities
844###############################################################################
845
846def SetScale(Data,refData):
847    Profile,Limits,Sample = Data
848    refProfile,refLimits,refSample = refData
849    x,y = Profile[:2]
850    rx,ry = refProfile[:2]
851    Beg = np.max([rx[0],x[0],Limits[1][0],refLimits[1][0]])
852    Fin = np.min([rx[-1],x[-1],Limits[1][1],refLimits[1][1]])
853    iBeg = np.searchsorted(x,Beg)
854    iFin = np.searchsorted(x,Fin)
855    sum = np.sum(y[iBeg:iFin])
856    refsum = np.sum(np.interp(x[iBeg:iFin],rx,ry,0,0))
857    Sample['Scale'][0] = refSample['Scale'][0]*refsum/sum
858   
859###############################################################################
860#### Size distribution
861###############################################################################
862
863def SizeDistribution(Profile,ProfDict,Limits,Substances,Sample,data):
864    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
865        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
866        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
867        'Unified disk':[UniDiskFF,UniDiskVol]}
868    Shape = data['Size']['Shape'][0]
869    Parms = data['Size']['Shape'][1:]
870    if data['Size']['logBins']:
871        Bins = np.logspace(np.log10(data['Size']['MinDiam']),np.log10(data['Size']['MaxDiam']),
872            data['Size']['Nbins']+1,True)/2.        #make radii
873    else:
874        Bins = np.linspace(data['Size']['MinDiam'],data['Size']['MaxDiam'],
875            data['Size']['Nbins']+1,True)/2.        #make radii
876    Dbins = np.diff(Bins)
877    Bins = Bins[:-1]+Dbins/2.
878    Contrast = Sample['Contrast'][1]
879    Scale = Sample['Scale'][0]
880    Sky = 10**data['Size']['MaxEnt']['Sky']
881    BinsBack = np.ones_like(Bins)*Sky*Scale/Contrast
882    Back = data['Back']
883    Q,Io,wt,Ic,Ib = Profile[:5]
884    Qmin = Limits[1][0]
885    Qmax = Limits[1][1]
886    wtFactor = ProfDict['wtFactor']
887    Ibeg = np.searchsorted(Q,Qmin)
888    Ifin = np.searchsorted(Q,Qmax)
889    BinMag = np.zeros_like(Bins)
890    Ic[:] = 0.
891    Gmat = G_matrix(Q[Ibeg:Ifin],Bins,Contrast,shapes[Shape][0],shapes[Shape][1],args=Parms)
892    if 'MaxEnt' == data['Size']['Method']:
893        chisq,BinMag,Ic[Ibeg:Ifin] = MaxEnt_SB(Scale*Io[Ibeg:Ifin]-Back[0],
894            Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),Gmat,BinsBack,
895            data['Size']['MaxEnt']['Niter'],report=True)
896    elif 'IPG' == data['Size']['Method']:
897        chisq,BinMag,Ic[Ibeg:Ifin] = IPG(Scale*Io[Ibeg:Ifin]-Back[0],Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),
898            Gmat,Bins,Dbins,data['Size']['IPG']['Niter'],Q[Ibeg:Ifin],approach=0.8,
899            Power=data['Size']['IPG']['Power'],report=True)
900    Ib[:] = Back[0]
901    Ic[Ibeg:Ifin] += Back[0]
902    print ' Final chi^2: %.3f'%(chisq)
903    Vols = shapes[Shape][1](Bins,Parms)
904    data['Size']['Distribution'] = [Bins,Dbins,BinMag/(2.*Dbins)]
905       
906################################################################################
907#### Unified fit
908################################################################################
909
910def UnifiedFit(Profile,ProfDict,Limits,Substances,Sample,data):
911    print 'do unified fit'
912   
913def UnifiedFxn(Q,G,Rg,B,Rgcf,P,SQfxn,args=[]):
914    termA = G*np.exp(-((Q*Rg)**2)/3.)
915    termB = B*np.exp(-((Q*Rgcf)**2)/3.)
916    termC = (scsp.erf(Q*Rg/np.sqrt(6))**3/Q)**P
917    return SQfxn(Q,args)*termA+(termB*termC)
918   
919
920
921################################################################################
922#### Modelling
923################################################################################
924
925def ModelFit(Profile,ProfDict,Limits,Substances,Sample,data):
926    print 'do model fit'
927   
928def ModelFxn(Profile,ProfDict,Limits,Substances,Sample,sasdData):
929    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
930        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
931        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
932        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol]}
933    partData = sasdData['Particle']
934    rhoMat = Substances['Substances'][partData['Matrix']['Name']].get('XAnom density',0.0)
935    matFrac = partData['Matrix']['VolFrac'] #[value,flag]       
936    Scale = Sample['Scale']     #[value,flag]
937    Back = sasdData['Back']
938    Q,Io,wt,Ic,Ib = Profile[:5]
939    Qmin = Limits[1][0]
940    Qmax = Limits[1][1]
941    wtFactor = ProfDict['wtFactor']
942    Ibeg = np.searchsorted(Q,Qmin)
943    Ifin = np.searchsorted(Q,Qmax)
944    Ib[:] = Back[0]
945    Ic[:] = 0
946    Ic[Ibeg:Ifin] += Back[0]
947    Rbins = []
948    Dist = []
949    for level in partData['Levels']:
950        controls = level['Controls']
951        FFfxn = shapes[controls['FormFact']][0]
952        Volfxn = shapes[controls['FormFact']][1]
953        FFargs = []
954        for item in ['Aspect ratio','Length','Thickness','Diameter',]:
955            if item in controls['FFargs']: 
956                FFargs.append(controls['FFargs'][item][0])
957        distFxn = controls['DistType']
958        rho = Substances['Substances'][level['Controls']['Material']].get('XAnom density',0.0)
959        contrast = rho**2-rhoMat**2
960        parmDict = level[controls['DistType']]
961        rBins,dist = MakeDiamDist(controls['DistType'],controls['NumPoints'],controls['Cutoff'],parmDict)
962        Gmat = G_matrix(Q[Ibeg:Ifin],rBins,contrast,FFfxn,Volfxn,FFargs)
963        dist *= level[distFxn]['Volume'][0]
964        Ic[Ibeg:Ifin] += np.dot(Gmat.T,dist)*100.
965        Rbins.append(rBins)
966        Dist.append(dist)
967       
968    return Rbins,Dist
969             
970   
971def MakeDiamDist(DistName,nPoints,cutoff,parmDict):
972   
973    if 'LogNormal' in DistName:
974        distFxn = 'LogNormalDist'
975        cumeFxn = 'LogNormalCume'
976        pos = parmDict['MinSize'][0]
977        args = [parmDict['Mean'][0],parmDict['StdDev'][0]]
978        step = 4.*np.sqrt(np.exp(parmDict['StdDev'][0]**2)*(np.exp(parmDict['StdDev'][0]**2)-1.))
979        mode = parmDict['MinSize'][0]+parmDict['Mean'][0]/np.exp(parmDict['StdDev'][0]**2)
980        minX = 1. #pos
981        maxX = 1.e4 #np.min([mode+nPoints*step*40,1.e5])
982    elif 'Gauss' in DistName:
983        distFxn = 'GaussDist'
984        cumeFxn = 'GaussCume'
985        pos = parmDict['Mean'][0]
986        args = [parmDict['StdDev'][0]]
987        step = 0.02*parmDict['StdDev'][0]
988        mode = parmDict['Mean'][0]
989        minX = np.max([mode-4.*parmDict['StdDev'][0],1.])
990        maxX = np.min([mode+4.*parmDict['StdDev'][0],1.e5])
991    elif 'LSW' in DistName:
992        distFxn = 'LSWDist'
993        cumeFxn = 'LSWCume'
994        pos = parmDict['Mean'][0]
995        args = []
996        minX,maxX = [0.,2.*pos]
997    elif 'Schulz' in DistName:
998        distFxn = 'SchulzZimmDist'
999        cumeFxn = 'SchulzZimmCume'
1000        pos = parmDict['Mean'][0]
1001        args = [parmDict['StdDev'][0]]
1002        minX = np.max([1.,pos-4.*parmDict['StdDev'][0]])
1003        maxX = np.min([pos+4.*parmDict['StdDev'][0],1.e5])
1004    nP = 500
1005    X = np.linspace(minX,maxX,nP,True)
1006    cume = eval(cumeFxn+'(X,pos,args)')
1007    Imin = np.searchsorted(cume,cutoff)
1008    Imax = np.min([np.searchsorted(cume,1-cutoff),nP-1])
1009    startX =  X[Imin]
1010    endX = X[Imax]
1011#    Diam = np.linspace(startX,endX,3*nPoints,True)
1012    Diam = np.logspace(0.,5.,500,True)
1013    TCW = eval(cumeFxn+'(Diam,pos,args)')
1014    CumeTgts = np.linspace(cutoff,(1.-cutoff),nPoints,True)
1015    rBins = np.interp(CumeTgts,TCW,Diam,0,0)
1016    return rBins,eval(distFxn+'(rBins,pos,args)')
1017
1018################################################################################
1019#### MaxEnt testing stuff
1020################################################################################
1021
1022def print_vec(text, a):
1023    '''print the contents of a vector to the console'''
1024    n = a.shape[0]
1025    print "%s[ = (" % text,
1026    for i in range(n):
1027        s = " %g, " % a[i]
1028        print s,
1029    print ")"
1030
1031def print_arr(text, a):
1032    '''print the contents of an array to the console'''
1033    n, m = a.shape
1034    print "%s[][] = (" % text
1035    for i in range(n):
1036        print " (",
1037        for j in range(m):
1038            print " %g, " % a[i][j],
1039        print "),"
1040    print ")"
1041
1042def test_MaxEnt_SB(report=True):
1043    def readTextData(filename):
1044        '''return q, I, dI from a 3-column text file'''
1045        if not os.path.exists(filename):
1046            raise Exception("file not found: " + filename)
1047        buf = [line.split() for line in open(filename, 'r').readlines()]
1048        M = len(buf)
1049        buf = zip(*buf)         # transpose rows and columns
1050        q  = np.array(buf[0], dtype=np.float64)
1051        I  = np.array(buf[1], dtype=np.float64)
1052        dI = np.array(buf[2], dtype=np.float64)
1053        return q, I, dI
1054    print "MaxEnt_SB: "
1055    test_data_file = os.path.join( 'testinp', 'test.sas')
1056    rhosq = 100     # scattering contrast, 10^20 1/cm^-4
1057    bkg   = 0.1     #   I = I - bkg
1058    dMin, dMax, nRadii = 25, 9000, 40
1059    defaultDistLevel = 1.0e-6
1060    IterMax = 40
1061    errFac = 1.05
1062   
1063    r    = np.logspace(math.log10(dMin), math.log10(dMax), nRadii)/2
1064    dr   = r * (r[1]/r[0] - 1)          # step size
1065    f_dr = np.ndarray((nRadii)) * 0  # volume fraction histogram
1066    b    = np.ndarray((nRadii)) * 0 + defaultDistLevel  # MaxEnt "sky background"
1067   
1068    qVec, I, dI = readTextData(test_data_file)
1069    G = G_matrix(qVec,r,rhosq,SphereFF,SphereVol,args=())
1070   
1071    chisq,f_dr,Ic = MaxEnt_SB(I - bkg, dI*errFac, b, IterMax, G, report=report)
1072    if f_dr is None:
1073        print "no solution"
1074        return
1075   
1076    print "solution reached"
1077    for a,b,c in zip(r.tolist(), dr.tolist(), f_dr.tolist()):
1078        print '%10.4f %10.4f %12.4g'%(a,b,c)
1079
1080def tests():
1081    test_MaxEnt_SB(report=True)
1082
1083if __name__ == '__main__':
1084    tests()
1085
Note: See TracBrowser for help on using the repository browser.