source: trunk/GSASIIsasd.py @ 1263

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

comment out the LEAVE_WINDOW in ValidatedTextCtrl? - leads to weird behavior when just moving mouse around.
Add SASD Model Add facility for Particle Modeling
Add Neutron CW SASD import & change default file extensions for SASD.
Begin Unified & Model fit routines - they don't do anything yet

File size: 32.1 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##### SB-MaxEnt
316################################################################################
317
318def G_matrix(q,r,contrast,FFfxn,Volfxn,args=()):
319    '''Calculates the response matrix :math:`G(Q,r)`
320   
321    :param float q: :math:`Q`
322    :param float r: :math:`r`
323    :param float contrast: :math:`|\\Delta\\rho|^2`, the scattering contrast
324    :param function FFfxn: form factor function FF(q,r,args)
325    :param function Volfxn: volume function Vol(r,args)
326    :returns float: G(Q,r)
327    '''
328    FF = FFfxn(q,r,args)
329    Vol = Volfxn(r,args)
330    return 1.e-4*(contrast*Vol*FF**2).T     #10^-20 vs 10^-24
331   
332'''
333sbmaxent
334
335Entropy maximization routine as described in the article
336J Skilling and RK Bryan; MNRAS 211 (1984) 111 - 124.
337("MNRAS": "Monthly Notices of the Royal Astronomical Society")
338
339:license: Copyright (c) 2013, UChicago Argonne, LLC
340:license: This file is distributed subject to a Software License Agreement found
341     in the file LICENSE that is included with this distribution.
342
343References:
344
3451. J Skilling and RK Bryan; MON NOT R ASTR SOC 211 (1984) 111 - 124.
3462. JA Potton, GJ Daniell, and BD Rainford; Proc. Workshop
347   Neutron Scattering Data Analysis, Rutherford
348   Appleton Laboratory, UK, 1986; ed. MW Johnson,
349   IOP Conference Series 81 (1986) 81 - 86, Institute
350   of Physics, Bristol, UK.
3513. ID Culverwell and GP Clarke; Ibid. 87 - 96.
3524. JA Potton, GK Daniell, & BD Rainford,
353   J APPL CRYST 21 (1988) 663 - 668.
3545. JA Potton, GJ Daniell, & BD Rainford,
355   J APPL CRYST 21 (1988) 891 - 897.
356
357'''
358
359class MaxEntException(Exception): 
360    '''Any exception from this module'''
361    pass
362
363def MaxEnt_SB(datum, sigma, G, base, IterMax, image_to_data=None, data_to_image=None, report=False):
364    '''
365    do the complete Maximum Entropy algorithm of Skilling and Bryan
366   
367    :param float datum[]:
368    :param float sigma[]:
369    :param float[][] G: transformation matrix
370    :param float base[]:
371    :param int IterMax:
372    :param obj image_to_data: opus function (defaults to opus)
373    :param obj data_to_image: tropus function (defaults to tropus)
374   
375    :returns float[]: :math:`f(r) dr`
376    '''
377       
378    TEST_LIMIT        = 0.05                    # for convergence
379    CHI_SQR_LIMIT     = 0.01                    # maximum difference in ChiSqr for a solution
380    SEARCH_DIRECTIONS = 3                       # <10.  This code requires value = 3
381    RESET_STRAYS      = 1                       # was 0.001, correction of stray negative values
382    DISTANCE_LIMIT_FACTOR = 0.1                 # limitation on df to constrain runaways
383   
384    MAX_MOVE_LOOPS    = 500                     # for no solution in routine: move,
385    MOVE_PASSES       = 0.001                   # convergence test in routine: move
386
387    def tropus (data, G):
388        '''
389        tropus: transform data-space -> solution-space:  [G] * data
390       
391        default definition, caller can use this definition or provide an alternative
392       
393        :param float[M] data: observations, ndarray of shape (M)
394        :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
395        :returns float[N]: calculated image, ndarray of shape (N)
396        '''
397        return G.dot(data)
398
399    def opus (image, G):
400        '''
401        opus: transform solution-space -> data-space:  [G]^tr * image
402       
403        default definition, caller can use this definition or provide an alternative
404       
405        :param float[N] image: solution, ndarray of shape (N)
406        :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
407        :returns float[M]: calculated data, ndarray of shape (M)
408        '''
409        return np.dot(G.T,image)    #G.transpose().dot(image)
410
411    def Dist(s2, beta):
412        '''measure the distance of this possible solution'''
413        w = 0
414        n = beta.shape[0]
415        for k in range(n):
416            z = -sum(s2[k] * beta)
417            w += beta[k] * z
418        return w
419   
420    def ChiNow(ax, c1, c2, s1, s2):
421        '''
422        ChiNow
423       
424        :returns tuple: (ChiNow computation of ``w``, beta)
425        '''
426       
427        bx = 1 - ax
428        a =   bx * c2  -  ax * s2
429        b = -(bx * c1  -  ax * s1)
430   
431        beta = ChoSol(a, b)
432        w = 1.0
433        for k in range(SEARCH_DIRECTIONS):
434            w += beta[k] * (c1[k] + 0.5*sum(c2[k] * beta))
435        return w, beta
436   
437    def ChoSol(a, b):
438        '''
439        ChoSol: ? chop the solution vectors ?
440       
441        :returns: new vector beta
442        '''
443        n = b.shape[0]
444        fl = np.zeros((n,n))
445        bl = np.zeros_like(b)
446       
447        #print_arr("ChoSol: a", a)
448        #print_vec("ChoSol: b", b)
449   
450        if (a[0][0] <= 0):
451            msg = "ChoSol: a[0][0] = " 
452            msg += str(a[0][0])
453            msg += '  Value must be positive'
454            raise MaxEntException(msg)
455   
456        # first, compute fl from a
457        # note fl is a lower triangular matrix
458        fl[0][0] = math.sqrt (a[0][0])
459        for i in (1, 2):
460            fl[i][0] = a[i][0] / fl[0][0]
461            for j in range(1, i+1):
462                z = 0.0
463                for k in range(j):
464                    z += fl[i][k] * fl[j][k]
465                    #print "ChoSol: %d %d %d  z = %lg" % ( i, j, k, z)
466                z = a[i][j] - z
467                if j == i:
468                    y = math.sqrt(z)
469                else:
470                    y = z / fl[j][j]
471                fl[i][j] = y
472        #print_arr("ChoSol: fl", fl)
473   
474        # next, compute bl from fl and b
475        bl[0] = b[0] / fl[0][0]
476        for i in (1, 2):
477            z = 0.0
478            for k in range(i):
479                z += fl[i][k] * bl[k]
480                #print "\t", i, k, z
481            bl[i] = (b[i] - z) / fl[i][i]
482        #print_vec("ChoSol: bl", bl)
483   
484        # last, compute beta from bl and fl
485        beta = np.empty((n))
486        beta[-1] = bl[-1] / fl[-1][-1]
487        for i in (1, 0):
488            z = 0.0
489            for k in range(i+1, n):
490                z += fl[k][i] * beta[k]
491                #print "\t\t", i, k, 'z=', z
492            beta[i] = (bl[i] - z) / fl[i][i]
493        #print_vec("ChoSol: beta", beta)
494   
495        return beta
496
497    def MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2):
498        '''
499        move beta one step closer towards the solution
500        '''
501        a_lower, a_upper = 0., 1.          # bracket  "a"
502        cmin, beta = ChiNow (a_lower, c1, c2, s1, s2)
503        #print "MaxEntMove: cmin = %g" % cmin
504        if cmin*chisq > chizer:
505            ctarg = (1.0 + cmin)/2
506        else:
507            ctarg = chizer/chisq
508        f_lower = cmin - ctarg
509        c_upper, beta = ChiNow (a_upper, c1, c2, s1, s2)
510        f_upper = c_upper - ctarg
511   
512        fx = 2*MOVE_PASSES      # just to start off
513        loop = 1
514        while abs(fx) >= MOVE_PASSES and loop <= MAX_MOVE_LOOPS:
515            a_new = (a_lower + a_upper) * 0.5           # search by bisection
516            c_new, beta = ChiNow (a_new, c1, c2, s1, s2)
517            fx = c_new - ctarg
518            # tighten the search range for the next pass
519            if f_lower*fx > 0:
520                a_lower, f_lower = a_new, fx
521            if f_upper*fx > 0:
522                a_upper, f_upper = a_new, fx
523            loop += 1
524   
525        if abs(fx) >= MOVE_PASSES or loop > MAX_MOVE_LOOPS:
526            msg = "MaxEntMove: Loop counter = " 
527            msg += str(MAX_MOVE_LOOPS)
528            msg += '  No convergence in alpha chop'
529            raise MaxEntException(msg)
530   
531        w = Dist (s2, beta);
532        m = SEARCH_DIRECTIONS
533        if (w > DISTANCE_LIMIT_FACTOR*fSum/blank):        # invoke the distance penalty, SB eq. 17
534            for k in range(m):
535                beta[k] *= math.sqrt (fSum/(blank*w))
536        chtarg = ctarg * chisq
537        return w, chtarg, loop, a_new, fx, beta
538       
539#MaxEnt_SB starts here   
540
541    if image_to_data == None:
542        image_to_data = opus
543    if data_to_image == None:
544        data_to_image = tropus
545    n   = len(base)
546    npt = len(datum)
547
548    # Note that the order of subscripts for
549    # "xi" and "eta" has been reversed from
550    # the convention used in the FORTRAN version
551    # to enable parts of them to be passed as
552    # as vectors to "image_to_data" and "data_to_image".
553    xi      = np.zeros((SEARCH_DIRECTIONS, n))
554    eta     = np.zeros((SEARCH_DIRECTIONS, npt))
555    beta    = np.zeros((SEARCH_DIRECTIONS))
556    # s1      = np.zeros((SEARCH_DIRECTIONS))
557    # c1      = np.zeros((SEARCH_DIRECTIONS))
558    s2      = np.zeros((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
559    c2      = np.zeros((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
560
561    # TODO: replace blank (scalar) with base (vector)
562    blank = sum(base) / len(base)   # use the average value of base
563
564    chizer, chtarg = npt*1.0, npt*1.0
565    f = base * 1.0                              # starting distribution is base
566
567    fSum  = sum(f)                              # find the sum of the f-vector
568    z = (datum - image_to_data (f, G)) / sigma  # standardized residuals, SB eq. 3
569    chisq = sum(z*z)                            # Chi^2, SB eq. 4
570
571    for iter in range(IterMax):
572        ox = -2 * z / sigma                        # gradient of Chi^2
573
574        cgrad = data_to_image (ox, G)              # cgrad[i] = del(C)/del(f[i]), SB eq. 8
575        sgrad = -np.log(f/base) / (blank*math.exp (1.0))  # sgrad[i] = del(S)/del(f[i])
576        snorm = math.sqrt(sum(f * sgrad*sgrad))    # entropy term, SB eq. 22
577        cnorm = math.sqrt(sum(f * cgrad*cgrad))    # ChiSqr term, SB eq. 22
578        tnorm = sum(f * sgrad * cgrad)             # norm for gradient term TEST
579
580        a = 1.0
581        b = 1.0 / cnorm
582        if iter == 0:
583            test = 0.0     # mismatch between entropy and ChiSquared gradients
584        else:
585            test = math.sqrt( ( 1.0 - tnorm/(snorm*cnorm) )/2 ) # SB eq. 37?
586            a = 0.5 / (snorm * test)
587            b *= 0.5 / test
588        xi[0] = f * cgrad / cnorm
589        xi[1] = f * (a * sgrad - b * cgrad)
590
591        eta[0] = image_to_data (xi[0], G);          # image --> data
592        eta[1] = image_to_data (xi[1], G);          # image --> data
593        ox = eta[1] / (sigma * sigma)
594        xi[2] = data_to_image (ox, G);              # data --> image
595        a = 1.0 / math.sqrt(sum(f * xi[2]*xi[2]))
596        xi[2] = f * xi[2] * a
597        eta[2] = image_to_data (xi[2], G)           # image --> data
598       
599#         print_arr("MaxEnt: eta.transpose()", eta.transpose())
600#         print_arr("MaxEnt: xi.transpose()", xi.transpose())
601
602        # prepare the search directions for the conjugate gradient technique
603        c1 = xi.dot(cgrad) / chisq                          # C_mu, SB eq. 24
604        s1 = xi.dot(sgrad)                          # S_mu, SB eq. 24
605#         print_vec("MaxEnt: c1", c1)
606#         print_vec("MaxEnt: s1", s1)
607
608        for k in range(SEARCH_DIRECTIONS):
609            for l in range(k+1):
610                c2[k][l] = 2 * sum(eta[k] * eta[l] / sigma/sigma) / chisq
611                s2[k][l] = -sum(xi[k] * xi[l] / f) / blank
612#         print_arr("MaxEnt: c2", c2)
613#         print_arr("MaxEnt: s2", s2)
614
615        # reflect across the body diagonal
616        for k, l in ((0,1), (0,2), (1,2)):
617            c2[k][l] = c2[l][k]                     #  M_(mu,nu)
618            s2[k][l] = s2[l][k]                     #  g_(mu,nu)
619 
620        beta[0] = -0.5 * c1[0] / c2[0][0]
621        beta[1] = 0.0
622        beta[2] = 0.0
623        if (iter > 0):
624            w, chtarg, loop, a_new, fx, beta = MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2)
625
626        f_old = f.copy()    # preserve the last image
627        f += xi.transpose().dot(beta)   # move the image towards the solution, SB eq. 25
628       
629        # As mentioned at the top of p.119,
630        # need to protect against stray negative values.
631        # In this case, set them to RESET_STRAYS * base[i]
632        #f = f.clip(RESET_STRAYS * blank, f.max())
633        for i in range(n):
634            if f[i] <= 0.0:
635                f[i] = RESET_STRAYS * base[i]
636        df = f - f_old
637        fSum = sum(f)
638        fChange = sum(df)
639
640        # calculate the normalized entropy
641        S = sum((f/fSum) * np.log(f/fSum))      # normalized entropy, S&B eq. 1
642        z = (datum - image_to_data (f, G)) / sigma  # standardized residuals
643        chisq = sum(z*z)                            # report this ChiSq
644
645        if report:
646            print " MaxEnt trial/max: %3d/%3d" % ((iter+1), IterMax)
647            print " Residual: %5.2lf%% Entropy: %8lg" % (100*test, S)
648            if iter > 0:
649                value = 100*( math.sqrt(chisq/chtarg)-1)
650            else:
651                value = 0
652      #      print " %12.5lg %10.4lf" % ( math.sqrt(chtarg/npt), value )
653            print " Function sum: %.6lg Change from last: %.2lf%%\n" % (fSum,100*fChange/fSum)
654
655        # See if we have finished our task.
656        # do the hardest test first
657        if (abs(chisq/chizer-1.0) < CHI_SQR_LIMIT) and  (test < TEST_LIMIT):
658            print ' Convergence achieved.'
659            return chisq,f,image_to_data(f, G)     # solution FOUND returns here
660    print ' No convergence! Try increasing Error multiplier.'
661    return chisq,f,image_to_data(f, G)       # no solution after IterMax iterations
662
663   
664###############################################################################
665#### IPG/TNNLS Routines
666###############################################################################
667
668def IPG(datum,sigma,G,Bins,Dbins,IterMax,Qvec=[],approach=0.8,Power=-1,report=False):
669    ''' An implementation of the Interior-Point Gradient method of
670    Michael Merritt & Yin Zhang, Technical Report TR04-08, Dept. of Comp. and
671    Appl. Math., Rice Univ., Houston, Texas 77005, U.S.A. found on the web at
672    http://www.caam.rice.edu/caam/trs/2004/TR04-08.pdf
673    Problem addressed: Total Non-Negative Least Squares (TNNLS)
674    :param float datum[]:
675    :param float sigma[]:
676    :param float[][] G: transformation matrix
677    :param int IterMax:
678    :param float Qvec: data positions for Power = 0-4
679    :param float approach: 0.8 default fitting parameter
680    :param int Power: 0-4 for Q^Power weighting, -1 to use input sigma
681   
682    '''
683    if Power < 0: 
684        GmatE = G/sigma[:np.newaxis]
685        IntE = datum/sigma
686        pwr = 0
687        QvecP = np.ones_like(datum)
688    else:
689        GmatE = G[:]
690        IntE = datum[:]
691        pwr = Power
692        QvecP = Qvec**pwr
693    Amat = GmatE*QvecP[:np.newaxis]
694    AAmat = np.inner(Amat,Amat)
695    Bvec = IntE*QvecP
696    Xw = np.ones_like(Bins)*1.e-32
697    calc = np.dot(G.T,Xw)
698    nIter = 0
699    err = 10.
700    while (nIter<IterMax) and (err > 1.):
701        #Step 1 in M&Z paper:
702        Qk = np.inner(AAmat,Xw)-np.inner(Amat,Bvec)
703        Dk = Xw/np.inner(AAmat,Xw)
704        Pk = -Dk*Qk
705        #Step 2 in M&Z paper:
706        alpSt = -np.inner(Pk,Qk)/np.inner(Pk,np.inner(AAmat,Pk))
707        alpWv = -Xw/Pk
708        AkSt = [np.where(Pk[i]<0,np.min((approach*alpWv[i],alpSt)),alpSt) for i in range(Pk.shape[0])]
709        #Step 3 in M&Z paper:
710        shift = AkSt*Pk
711        Xw += shift
712        #done IPG; now results
713        nIter += 1
714        calc = np.dot(G.T,Xw)
715        chisq = np.sum(((datum-calc)/sigma)**2)
716        err = chisq/len(datum)
717        if report:
718            print ' Iteration: %d, chisq: %.3g, sum(shift^2): %.3g'%(nIter,chisq,np.sum(shift**2))
719    return chisq,Xw,calc
720
721###############################################################################
722#### SASD Utilities
723###############################################################################
724
725def SetScale(Data,refData):
726    Profile,Limits,Sample = Data
727    refProfile,refLimits,refSample = refData
728    x,y = Profile[:2]
729    rx,ry = refProfile[:2]
730    Beg = np.max([rx[0],x[0],Limits[1][0],refLimits[1][0]])
731    Fin = np.min([rx[-1],x[-1],Limits[1][1],refLimits[1][1]])
732    iBeg = np.searchsorted(x,Beg)
733    iFin = np.searchsorted(x,Fin)
734    sum = np.sum(y[iBeg:iFin])
735    refsum = np.sum(np.interp(x[iBeg:iFin],rx,ry,0,0))
736    Sample['Scale'][0] = refSample['Scale'][0]*refsum/sum
737   
738###############################################################################
739#### Size distribution
740###############################################################################
741
742def SizeDistribution(Profile,ProfDict,Limits,Substances,Sample,data):
743    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
744        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
745        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
746        'Unified disk':[UniDiskFF,UniDiskVol]}
747    Shape = data['Size']['Shape'][0]
748    Parms = data['Size']['Shape'][1:]
749    if data['Size']['logBins']:
750        Bins = np.logspace(np.log10(data['Size']['MinDiam']),np.log10(data['Size']['MaxDiam']),
751            data['Size']['Nbins']+1,True)/2.        #make radii
752    else:
753        Bins = np.linspace(data['Size']['MinDiam'],data['Size']['MaxDiam'],
754            data['Size']['Nbins']+1,True)/2.        #make radii
755    Dbins = np.diff(Bins)
756    Bins = Bins[:-1]+Dbins/2.
757    Contrast = Sample['Contrast'][1]
758    Scale = Sample['Scale'][0]
759    Sky = 10**data['Size']['MaxEnt']['Sky']
760    BinsBack = np.ones_like(Bins)*Sky*Scale/Contrast
761    Back = data['Back']
762    Q,Io,wt,Ic,Ib = Profile[:5]
763    Qmin = Limits[1][0]
764    Qmax = Limits[1][1]
765    wtFactor = ProfDict['wtFactor']
766    Ibeg = np.searchsorted(Q,Qmin)
767    Ifin = np.searchsorted(Q,Qmax)
768    BinMag = np.zeros_like(Bins)
769    Ic[:] = 0.
770    Gmat = G_matrix(Q[Ibeg:Ifin],Bins,Contrast,shapes[Shape][0],shapes[Shape][1],args=Parms)
771    if 'MaxEnt' == data['Size']['Method']:
772        chisq,BinMag,Ic[Ibeg:Ifin] = MaxEnt_SB(Scale*Io[Ibeg:Ifin]-Back[0],
773            Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),Gmat,BinsBack,
774            data['Size']['MaxEnt']['Niter'],report=True)
775    elif 'IPG' == data['Size']['Method']:
776        chisq,BinMag,Ic[Ibeg:Ifin] = IPG(Scale*Io[Ibeg:Ifin]-Back[0],Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),
777            Gmat,Bins,Dbins,data['Size']['IPG']['Niter'],Q[Ibeg:Ifin],approach=0.8,
778            Power=data['Size']['IPG']['Power'],report=True)
779    Ib[:] = Back[0]
780    Ic[Ibeg:Ifin] += Back[0]
781    print ' Final chi^2: %.3f'%(chisq)
782    Vols = shapes[Shape][1](Bins,Parms)
783    data['Size']['Distribution'] = [Bins,Dbins,BinMag/(2.*Dbins)]
784       
785################################################################################
786#### Unified fit
787################################################################################
788
789def UnifiedFit(Profile,ProfDict,Limits,Substances,Sample,data):
790    print 'do unified fit'
791   
792def UnifiedFxn(Q,G,Rg,B,Rgcf,P,SQfxn,args=[]):
793    termA = G*np.exp(-((Q*Rg)**2)/3.)
794    termB = B*np.exp(-((Q*Rgcf)**2)/3.)
795    termC = (scsp.erf(Q*Rg/np.sqrt(6))**3/Q)**P
796    return SQfxn(Q,args)*termA+(termB*termC)
797   
798
799
800################################################################################
801#### Modelling
802################################################################################
803
804def ModelFit(Profile,ProfDict,Limits,Substances,Sample,data):
805    print 'do model fit'
806           
807   
808################################################################################
809#### MaxEnt testing stuff
810################################################################################
811
812def print_vec(text, a):
813    '''print the contents of a vector to the console'''
814    n = a.shape[0]
815    print "%s[ = (" % text,
816    for i in range(n):
817        s = " %g, " % a[i]
818        print s,
819    print ")"
820
821def print_arr(text, a):
822    '''print the contents of an array to the console'''
823    n, m = a.shape
824    print "%s[][] = (" % text
825    for i in range(n):
826        print " (",
827        for j in range(m):
828            print " %g, " % a[i][j],
829        print "),"
830    print ")"
831
832def test_MaxEnt_SB(report=True):
833    def readTextData(filename):
834        '''return q, I, dI from a 3-column text file'''
835        if not os.path.exists(filename):
836            raise Exception("file not found: " + filename)
837        buf = [line.split() for line in open(filename, 'r').readlines()]
838        M = len(buf)
839        buf = zip(*buf)         # transpose rows and columns
840        q  = np.array(buf[0], dtype=np.float64)
841        I  = np.array(buf[1], dtype=np.float64)
842        dI = np.array(buf[2], dtype=np.float64)
843        return q, I, dI
844    print "MaxEnt_SB: "
845    test_data_file = os.path.join( 'testinp', 'test.sas')
846    rhosq = 100     # scattering contrast, 10^20 1/cm^-4
847    bkg   = 0.1     #   I = I - bkg
848    dMin, dMax, nRadii = 25, 9000, 40
849    defaultDistLevel = 1.0e-6
850    IterMax = 40
851    errFac = 1.05
852   
853    r    = np.logspace(math.log10(dMin), math.log10(dMax), nRadii)/2
854    dr   = r * (r[1]/r[0] - 1)          # step size
855    f_dr = np.ndarray((nRadii)) * 0  # volume fraction histogram
856    b    = np.ndarray((nRadii)) * 0 + defaultDistLevel  # MaxEnt "sky background"
857   
858    qVec, I, dI = readTextData(test_data_file)
859    G = G_matrix(qVec,r,rhosq,SphereFF,SphereVol,args=())
860   
861    chisq,f_dr,Ic = MaxEnt_SB(I - bkg, dI*errFac, b, IterMax, G, report=report)
862    if f_dr is None:
863        print "no solution"
864        return
865   
866    print "solution reached"
867    for a,b,c in zip(r.tolist(), dr.tolist(), f_dr.tolist()):
868        print '%10.4f %10.4f %12.4g'%(a,b,c)
869
870def tests():
871    test_MaxEnt_SB(report=True)
872
873if __name__ == '__main__':
874    tests()
875
Note: See TracBrowser for help on using the repository browser.