source: trunk/GSASIIsasd.py @ 1271

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

remove a bit of dead code

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