source: trunk/GSASIIsasd.py @ 1266

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

change some formats in Image Controls
enable insertion of additional points in polygon & frame masks
add code for particle distributions

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