source: trunk/GSASIIsasd.py @ 1246

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

add PlotSASDSizeDist to auto plots for SASD/Model
fix shapes form factors & volumes to work with numpy arrays - all seem to work

File size: 26.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 sys
16import math
17import time
18
19import numpy as np
20import scipy as sp
21import numpy.linalg as nl
22from numpy.fft import ifft, fft, fftshift
23import scipy.special as scsp
24import scipy.interpolate as si
25import scipy.stats as st
26import scipy.optimize as so
27
28import GSASIIpath
29GSASIIpath.SetVersionNumber("$Revision: 1186 $")
30import GSASIIlattice as G2lat
31import GSASIIspc as G2spc
32import GSASIIElem as G2elem
33import GSASIIgrid as G2gd
34import GSASIIIO as G2IO
35import GSASIImath as G2mth
36import pypowder as pyd
37
38# trig functions in degrees
39sind = lambda x: math.sin(x*math.pi/180.)
40asind = lambda x: 180.*math.asin(x)/math.pi
41tand = lambda x: math.tan(x*math.pi/180.)
42atand = lambda x: 180.*math.atan(x)/math.pi
43atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
44cosd = lambda x: math.cos(x*math.pi/180.)
45acosd = lambda x: 180.*math.acos(x)/math.pi
46rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
47#numpy versions
48npsind = lambda x: np.sin(x*np.pi/180.)
49npasind = lambda x: 180.*np.arcsin(x)/math.pi
50npcosd = lambda x: np.cos(x*math.pi/180.)
51npacosd = lambda x: 180.*np.arccos(x)/math.pi
52nptand = lambda x: np.tan(x*math.pi/180.)
53npatand = lambda x: 180.*np.arctan(x)/np.pi
54npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
55npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave
56npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
57   
58###############################################################################
59#### Particle form factors
60###############################################################################
61
62def SphereFF(Q,R,args=()):
63    ''' Compute hard sphere form factor - can use numpy arrays
64    param float:Q Q value array (usually in A-1)
65    param float:R sphere radius (Usually in A - must match Q-1 units)
66    returns float: form factors as array as needed
67    '''
68    QR = Q[:,np.newaxis]*R
69    return (3./(QR**3))*(np.sin(QR)-(QR*np.cos(QR)))
70   
71def SpheroidFF(Q,R,args):
72    ''' Compute form factor of cylindrically symmetric ellipsoid (spheroid)
73    - can use numpy arrays for R & AR; will return corresponding numpy array
74    param float:Q Q value array (usually in A-1)
75    param float R: radius along 2 axes of spheroid
76    param float AR: aspect ratio so 3rd axis = R*AR
77    returns float: form factors as array as needed
78    '''
79    NP = 50
80    AR = args[0]
81    if 0.99 < AR < 1.01:
82        return SphereFF(Q,R,0)
83    else:
84        cth = np.linspace(0,1.,NP)
85        Rct = R[:,np.newaxis]*np.sqrt(1.+(AR**2-1.)*cth**2)
86        return np.sqrt(np.sum(SphereFF(Q[:,np.newaxis],Rct,0)**2,axis=2)/NP)
87           
88def CylinderFF(Q,R,args):
89    ''' Compute form factor for cylinders - can use numpy arrays
90    param float: Q Q value array (A-1)
91    param float: R cylinder radius (A)
92    param float: L cylinder length (A)
93    returns float: form factor
94    '''
95    L = args[0]
96    NP = 200
97    alp = np.linspace(0,np.pi/2.,NP)
98    LBessArg = 0.5*L*(Q[:,np.newaxis]*np.cos(alp))
99    LBess = np.where(LBessArg<1.e-6,1.,np.sin(LBessArg)/LBessArg)
100    QR = Q[:,np.newaxis]*R
101    SBessArg = QR[:,:,np.newaxis]*np.sin(alp)
102    SBess = np.where(SBessArg<1.e-6,0.5,scsp.jv(1,SBessArg)/SBessArg)
103    LSBess = LBess[:,np.newaxis,:]*SBess
104    return np.sqrt(2.*np.pi*np.sum(np.sin(alp)*LSBess**2,axis=2)/NP)
105   
106def CylinderDFF(Q,L,args):
107    ''' Compute form factor for cylinders - can use numpy arrays
108    param float: Q Q value array (A-1)
109    param float: L cylinder half length (A)
110    param float: R cylinder diameter (A)
111    returns float: form factor
112    '''
113    R = args[0]
114    return CylinderFF(Q,R,2.*L)   
115   
116def CylinderARFF(Q,R,args): 
117    ''' Compute form factor for cylinders - can use numpy arrays
118    param float: Q Q value array (A-1)
119    param float: R cylinder radius (A)
120    param float: AR cylinder aspect ratio = L/D = L/2R
121    returns float: form factor
122    '''
123    AR = args[0]
124    return CylinderFF(Q,R,2.*R*AR)   
125   
126def UniSphereFF(Q,R,args=0):
127    Rg = np.sqrt(3./5.)*R
128    B = 1.62/(Rg**4)    #are we missing *np.pi? 1.62 = 6*(3/5)**2/(4/3) sense?
129    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg/np.sqrt(6)))**3
130    return np.sqrt(np.exp((-Q[:,np.newaxis]**2*Rg**2)/3.)+(B/QstV**4))
131   
132def UniRodFF(Q,R,args):
133    L = args[0]
134    Rg2 = np.sqrt(R**2/2+L**2/12)
135    B2 = np.pi/L
136    Rg1 = np.sqrt(3.)*R/2.
137    G1 = (2./3.)*R/L
138    B1 = 4.*(L+R)/(R**3*L**2)
139    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg2/np.sqrt(6)))**3
140    FF = np.exp(-Q[:,np.newaxis]**2*Rg2**2/3.)+(B2/QstV)*np.exp(-Rg1**2*Q[:,np.newaxis]**2/3.)
141    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg1/np.sqrt(6)))**3
142    FF += G1*np.exp(-Q[:,np.newaxis]**2*Rg1**2/3.)+(B1/QstV**4)
143    return np.sqrt(FF)
144   
145def UniRodARFF(Q,R,args):
146    AR = args[0]
147    return UniRodFF(Q,R,[AR*R,])
148   
149def UniDiskFF(Q,R,args):
150    T = args[0]
151    Rg2 = np.sqrt(R**2/2.+T**2/12.)
152    B2 = 2./R**2
153    Rg1 = np.sqrt(3.)*T/2.
154    RgC2 = 1.1*Rg1
155    G1 = (2./3.)*(T/R)**2
156    B1 = 4.*(T+R)/(R**3*T**2)
157    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg2/np.sqrt(6)))**3
158    FF = np.exp(-Q[:,np.newaxis]**2*Rg2**2/3.)+(B2/QstV**2)*np.exp(-RgC2**2*Q[:,np.newaxis]**2/3.)
159    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg1/np.sqrt(6)))**3
160    FF += G1*np.exp(-Q[:,np.newaxis]**2*Rg1**2/3.)+(B1/QstV**4)
161    return np.sqrt(FF)
162   
163def UniTubeFF(Q,R,args):
164    L,T = args[:2]
165    Ri = R-T
166    DR2 = R**2-Ri**2
167    Vt = np.pi*DR2*L
168    Rg3 = np.sqrt(DR2/2.+L**2/12.)
169    B1 = 4.*np.pi**2*(DR2+L*(R+Ri))/Vt**2
170    B2 = np.pi**2*T/Vt
171    B3 = np.pi/L
172    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg3/np.sqrt(6)))**3
173    FF = np.exp(-Q[:,np.newaxis]**2*Rg3**2/3.)+(B3/QstV)*np.exp(-Q[:,np.newaxis]**2*R**2/3.)
174    QstV = Q/(scsp.erf(Q[:,np.newaxis]*R/np.sqrt(6)))**3
175    FF += (B2/QstV**2)*np.exp(-Q[:,np.newaxis]**2*T**2/3.)
176    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*T/np.sqrt(6)))**3
177    FF += B1/QstV**4
178    return np.sqrt(FF)
179
180###############################################################################
181#### Particle volumes
182###############################################################################
183
184def SphereVol(R,args=()):
185    ''' Compute volume of sphere
186    - numpy array friendly
187    param float:R sphere radius
188    returns float: volume
189    '''
190    return (4./3.)*np.pi*R**3
191
192def SpheroidVol(R,args):
193    ''' Compute volume of cylindrically symmetric ellipsoid (spheroid)
194    - numpy array friendly
195    param float R: radius along 2 axes of spheroid
196    param float AR: aspect ratio so radius of 3rd axis = R*AR
197    returns float: volume
198    '''
199    AR = args[0]
200    return AR*SphereVol(R)
201   
202def CylinderVol(R,args):
203    ''' Compute cylinder volume for radius & length
204    - numpy array friendly
205    param float: R diameter (A)
206    param float: L length (A)
207    returns float:volume (A^3)
208    '''
209    L = args[0]
210    return np.pi*L*R**2
211   
212def CylinderDVol(L,args):
213    ''' Compute cylinder volume for length & diameter
214    - numpy array friendly
215    param float: L half length (A)
216    param float: D diameter (A)
217    returns float:volume (A^3)
218    '''
219    D = args[0]
220    return CylinderVol(D/2.,[2.*L,])
221   
222def CylinderARVol(R,args):
223    ''' Compute cylinder volume for radius & aspect ratio = L/D
224    - numpy array friendly
225    param float: R radius (A)
226    param float: AR=L/D=L/2R aspect ratio
227    returns float:volume
228    '''
229    AR = args[0]
230    return CylinderVol(R,[2.*R*AR,])
231   
232def UniSphereVol(R,args=()):
233    ''' Compute volume of sphere
234    - numpy array friendly
235    param float:R sphere radius
236    returns float: volume
237    '''
238    return SphereVol(R)
239   
240def UniRodVol(R,args):
241    ''' Compute cylinder volume for radius & length
242    - numpy array friendly
243    param float: R diameter (A)
244    param float: L length (A)
245    returns float:volume (A^3)
246    '''
247    L = args[0]
248    return CylinderVol(R,[L,])
249   
250def UniRodARVol(R,args):
251    AR = args[0]
252    return CylinderARVol(R,[AR,])
253   
254def UniDiskVol(R,args):
255    T = args[0]
256    return CylinderVol(R,[T,])
257   
258def UniTubeVol(R,args):
259    ''' Compute tube volume for radius, length & wall thickness
260    - numpy array friendly
261    param float: R diameter (A)
262    param float: L length (A)
263    param float: T tube wall thickness (A)
264    returns float: volume (A^3) of tube wall
265    '''
266    L,T = arg[:2]
267    return CylinderVol(R,[L,])-CylinderVol(R-T,[L,])
268         
269################################################################################
270##### SB-MaxEnt
271################################################################################
272
273def G_matrix(q,r,contrast,FFfxn,Volfxn,args=()):
274    '''Calculates the response matrix :math:`G(Q,r)`
275   
276    :param float q: :math:`Q`
277    :param float r: :math:`r`
278    :param float contrast: :math:`|\\Delta\\rho|^2`, the scattering contrast
279    :param function FFfxn: form factor function FF(q,r,args)
280    :param function Volfxn: volume function Vol(r,args)
281    :returns float: G(Q,r)
282    '''
283    FF = FFfxn(q,r,args)
284    Vol = Volfxn(r,args)
285    return 1.e-4*(contrast*Vol*FF**2)     #10^-20 vs 10^-24
286   
287'''
288sbmaxent
289
290Entropy maximization routine as described in the article
291J Skilling and RK Bryan; MNRAS 211 (1984) 111 - 124.
292("MNRAS": "Monthly Notices of the Royal Astronomical Society")
293
294:license: Copyright (c) 2013, UChicago Argonne, LLC
295:license: This file is distributed subject to a Software License Agreement found
296     in the file LICENSE that is included with this distribution.
297
298References:
299
3001. J Skilling and RK Bryan; MON NOT R ASTR SOC 211 (1984) 111 - 124.
3012. JA Potton, GJ Daniell, and BD Rainford; Proc. Workshop
302   Neutron Scattering Data Analysis, Rutherford
303   Appleton Laboratory, UK, 1986; ed. MW Johnson,
304   IOP Conference Series 81 (1986) 81 - 86, Institute
305   of Physics, Bristol, UK.
3063. ID Culverwell and GP Clarke; Ibid. 87 - 96.
3074. JA Potton, GK Daniell, & BD Rainford,
308   J APPL CRYST 21 (1988) 663 - 668.
3095. JA Potton, GJ Daniell, & BD Rainford,
310   J APPL CRYST 21 (1988) 891 - 897.
311
312'''
313
314import os
315import sys
316import math
317import numpy
318
319class MaxEntException(Exception): 
320    '''Any exception from this module'''
321    pass
322
323def MaxEnt_SB(datum, sigma, base, IterMax, G, image_to_data=None, data_to_image=None, report=False):
324    '''
325    do the complete Maximum Entropy algorithm of Skilling and Bryan
326   
327    :param float datum[]:
328    :param float sigma[]:
329    :param float base[]:
330    :param int IterMax:
331    :param float[][] G: transformation matrix
332    :param obj image_to_data: opus function (defaults to opus)
333    :param obj data_to_image: tropus function (defaults to tropus)
334   
335    :returns float[]: :math:`f(r) dr`
336    '''
337       
338    TEST_LIMIT        = 0.05                    # for convergence
339    CHI_SQR_LIMIT     = 0.01                    # maximum difference in ChiSqr for a solution
340    SEARCH_DIRECTIONS = 3                       # <10.  This code requires value = 3
341    RESET_STRAYS      = 1                       # was 0.001, correction of stray negative values
342    DISTANCE_LIMIT_FACTOR = 0.1                 # limitation on df to constrain runaways
343   
344    MAX_MOVE_LOOPS    = 500                     # for no solution in routine: move,
345    MOVE_PASSES       = 0.001                   # convergence test in routine: move
346
347    def opus (data, G):
348        '''
349        opus: transform data-space -> solution-space:  [G] * data
350       
351        default definition, caller can use this definition or provide an alternative
352       
353        :param float[M] data: observations, ndarray of shape (M)
354        :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
355        :returns float[N]: calculated image, ndarray of shape (N)
356        '''
357        return G.dot(data)
358
359    def tropus (image, G):
360        '''
361        tropus: transform solution-space -> data-space:  [G]^tr * image
362       
363        default definition, caller can use this definition or provide an alternative
364       
365        :param float[N] image: solution, ndarray of shape (N)
366        :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
367        :returns float[M]: calculated data, ndarray of shape (M)
368        '''
369        return G.transpose().dot(image)
370
371    def Dist(s2, beta):
372        '''measure the distance of this possible solution'''
373        w = 0
374        n = beta.shape[0]
375        for k in range(n):
376            z = -sum(s2[k] * beta)
377            w += beta[k] * z
378        return w
379   
380    def ChiNow(ax, c1, c2, s1, s2):
381        '''
382        ChiNow
383       
384        :returns tuple: (ChiNow computation of ``w``, beta)
385        '''
386       
387        bx = 1 - ax
388        a =   bx * c2  -  ax * s2
389        b = -(bx * c1  -  ax * s1)
390   
391        beta = ChoSol(a, b)
392        w = 1.0
393        for k in range(SEARCH_DIRECTIONS):
394            w += beta[k] * (c1[k] + 0.5*sum(c2[k] * beta))
395        return w, beta
396   
397    def ChoSol(a, b):
398        '''
399        ChoSol: ? chop the solution vectors ?
400       
401        :returns: new vector beta
402        '''
403        n = b.shape[0]
404        fl = numpy.ndarray((n, n))*0
405        bl = numpy.ndarray((n))*0
406       
407        #print_arr("ChoSol: a", a)
408        #print_vec("ChoSol: b", b)
409   
410        if (a[0][0] <= 0):
411            msg = "ChoSol: a[0][0] = " 
412            msg += str(a[0][0])
413            msg += '  Value must be positive'
414            raise MaxEntException(msg)
415   
416        # first, compute fl from a
417        # note fl is a lower triangular matrix
418        fl[0][0] = math.sqrt (a[0][0])
419        for i in (1, 2):
420            fl[i][0] = a[i][0] / fl[0][0]
421            for j in range(1, i+1):
422                z = 0.0
423                for k in range(j):
424                    z += fl[i][k] * fl[j][k]
425                    #print "ChoSol: %d %d %d  z = %lg" % ( i, j, k, z)
426                z = a[i][j] - z
427                if j == i:
428                    y = math.sqrt(z)
429                else:
430                    y = z / fl[j][j]
431                fl[i][j] = y
432        #print_arr("ChoSol: fl", fl)
433   
434        # next, compute bl from fl and b
435        bl[0] = b[0] / fl[0][0]
436        for i in (1, 2):
437            z = 0.0
438            for k in range(i):
439                z += fl[i][k] * bl[k]
440                #print "\t", i, k, z
441            bl[i] = (b[i] - z) / fl[i][i]
442        #print_vec("ChoSol: bl", bl)
443   
444        # last, compute beta from bl and fl
445        beta = numpy.ndarray((n))
446        beta[-1] = bl[-1] / fl[-1][-1]
447        for i in (1, 0):
448            z = 0.0
449            for k in range(i+1, n):
450                z += fl[k][i] * beta[k]
451                #print "\t\t", i, k, 'z=', z
452            beta[i] = (bl[i] - z) / fl[i][i]
453        #print_vec("ChoSol: beta", beta)
454   
455        return beta
456
457    def MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2):
458        '''
459        move beta one step closer towards the solution
460        '''
461        a_lower, a_upper = 0., 1.          # bracket  "a"
462        cmin, beta = ChiNow (a_lower, c1, c2, s1, s2)
463        #print "MaxEntMove: cmin = %g" % cmin
464        if cmin*chisq > chizer:
465            ctarg = (1.0 + cmin)/2
466        else:
467            ctarg = chizer/chisq
468        f_lower = cmin - ctarg
469        c_upper, beta = ChiNow (a_upper, c1, c2, s1, s2)
470        f_upper = c_upper - ctarg
471   
472        fx = 2*MOVE_PASSES      # just to start off
473        loop = 1
474        while abs(fx) >= MOVE_PASSES and loop <= MAX_MOVE_LOOPS:
475            a_new = (a_lower + a_upper) * 0.5           # search by bisection
476            c_new, beta = ChiNow (a_new, c1, c2, s1, s2)
477            fx = c_new - ctarg
478            # tighten the search range for the next pass
479            if f_lower*fx > 0:
480                a_lower, f_lower = a_new, fx
481            if f_upper*fx > 0:
482                a_upper, f_upper = a_new, fx
483            loop += 1
484   
485        if abs(fx) >= MOVE_PASSES or loop > MAX_MOVE_LOOPS:
486            msg = "MaxEntMove: Loop counter = " 
487            msg += str(MAX_MOVE_LOOPS)
488            msg += '  No convergence in alpha chop'
489            raise MaxEntException(msg)
490   
491        w = Dist (s2, beta);
492        m = SEARCH_DIRECTIONS
493        if (w > DISTANCE_LIMIT_FACTOR*fSum/blank):        # invoke the distance penalty, SB eq. 17
494            for k in range(m):
495                beta[k] *= math.sqrt (fSum/(blank*w))
496        chtarg = ctarg * chisq
497        return w, chtarg, loop, a_new, fx, beta
498       
499#MaxEnt_SB starts here   
500
501    if image_to_data == None:
502        image_to_data = opus
503    if data_to_image == None:
504        data_to_image = tropus
505    n   = len(base)
506    npt = len(datum)
507
508    # Note that the order of subscripts for
509    # "xi" and "eta" has been reversed from
510    # the convention used in the FORTRAN version
511    # to enable parts of them to be passed as
512    # as vectors to "image_to_data" and "data_to_image".
513    xi      = 0*numpy.ndarray((SEARCH_DIRECTIONS, n))
514    eta     = 0*numpy.ndarray((SEARCH_DIRECTIONS, npt))
515    beta    = 0*numpy.ndarray((SEARCH_DIRECTIONS))
516    # s1      = 0*numpy.ndarray((SEARCH_DIRECTIONS))
517    # c1      = 0*numpy.ndarray((SEARCH_DIRECTIONS))
518    s2      = 0*numpy.ndarray((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
519    c2      = 0*numpy.ndarray((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
520
521    # TODO: replace blank (scalar) with base (vector)
522    blank = sum(base) / len(base)   # use the average value of base
523
524    chizer, chtarg = npt*1.0, npt*1.0
525    f = base * 1.0                              # starting distribution is base
526
527    fSum  = sum(f)                              # find the sum of the f-vector
528    z = (datum - image_to_data (f, G)) / sigma  # standardized residuals, SB eq. 3
529    chisq = sum(z*z)                            # Chi^2, SB eq. 4
530
531    for iter in range(IterMax):
532        ox = -2 * z / sigma                        # gradient of Chi^2
533
534        cgrad = data_to_image (ox, G)              # cgrad[i] = del(C)/del(f[i]), SB eq. 8
535        sgrad = -numpy.log(f/base) / (blank*math.exp (1.0))  # sgrad[i] = del(S)/del(f[i])
536        snorm = math.sqrt(sum(f * sgrad*sgrad))    # entropy term, SB eq. 22
537        cnorm = math.sqrt(sum(f * cgrad*cgrad))    # ChiSqr term, SB eq. 22
538        tnorm = sum(f * sgrad * cgrad)             # norm for gradient term TEST
539
540        a = 1.0
541        b = 1.0 / cnorm
542        if iter == 0:
543            test = 0.0     # mismatch between entropy and ChiSquared gradients
544        else:
545            test = math.sqrt( ( 1.0 - tnorm/(snorm*cnorm) )/2 ) # SB eq. 37?
546            a = 0.5 / (snorm * test)
547            b *= 0.5 / test
548        xi[0] = f * cgrad / cnorm
549        xi[1] = f * (a * sgrad - b * cgrad)
550
551        eta[0] = image_to_data (xi[0], G);          # image --> data
552        eta[1] = image_to_data (xi[1], G);          # image --> data
553        ox = eta[1] / (sigma * sigma)
554        xi[2] = data_to_image (ox, G);              # data --> image
555        a = 1.0 / math.sqrt(sum(f * xi[2]*xi[2]))
556        xi[2] = f * xi[2] * a
557        eta[2] = image_to_data (xi[2], G)           # image --> data
558       
559#         print_arr("MaxEnt: eta.transpose()", eta.transpose())
560#         print_arr("MaxEnt: xi.transpose()", xi.transpose())
561
562        # prepare the search directions for the conjugate gradient technique
563        c1 = xi.dot(cgrad) / chisq                          # C_mu, SB eq. 24
564        s1 = xi.dot(sgrad)                          # S_mu, SB eq. 24
565#         print_vec("MaxEnt: c1", c1)
566#         print_vec("MaxEnt: s1", s1)
567
568        for k in range(SEARCH_DIRECTIONS):
569            for l in range(k+1):
570                c2[k][l] = 2 * sum(eta[k] * eta[l] / sigma/sigma) / chisq
571                s2[k][l] = -sum(xi[k] * xi[l] / f) / blank
572#         print_arr("MaxEnt: c2", c2)
573#         print_arr("MaxEnt: s2", s2)
574
575        # reflect across the body diagonal
576        for k, l in ((0,1), (0,2), (1,2)):
577            c2[k][l] = c2[l][k]                     #  M_(mu,nu)
578            s2[k][l] = s2[l][k]                     #  g_(mu,nu)
579 
580        beta[0] = -0.5 * c1[0] / c2[0][0]
581        beta[1] = 0.0
582        beta[2] = 0.0
583        if (iter > 0):
584            w, chtarg, loop, a_new, fx, beta = MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2)
585
586        f_old = f.copy()    # preserve the last image
587        f += xi.transpose().dot(beta)   # move the image towards the solution, SB eq. 25
588       
589        # As mentioned at the top of p.119,
590        # need to protect against stray negative values.
591        # In this case, set them to RESET_STRAYS * base[i]
592        #f = f.clip(RESET_STRAYS * blank, f.max())
593        for i in range(n):
594            if f[i] <= 0.0:
595                f[i] = RESET_STRAYS * base[i]
596        df = f - f_old
597        fSum = sum(f)
598        fChange = sum(df)
599
600        # calculate the normalized entropy
601        S = sum((f/fSum) * numpy.log(f/fSum))      # normalized entropy, S&B eq. 1
602        z = (datum - image_to_data (f, G)) / sigma  # standardized residuals
603        chisq = sum(z*z)                            # report this ChiSq
604
605        if report:
606            print " MaxEnt trial/max: %3d/%3d" % ((iter+1), IterMax)
607            print " Residual: %5.2lf%% Entropy: %8lg" % (100*test, S)
608            if iter > 0:
609                value = 100*( math.sqrt(chisq/chtarg)-1)
610            else:
611                value = 0
612      #      print " %12.5lg %10.4lf" % ( math.sqrt(chtarg/npt), value )
613            print " Function sum: %.6lg Change from last: %.2lf%%\n" % (fSum,100*fChange/fSum)
614
615        # See if we have finished our task.
616        # do the hardest test first
617        if (abs(chisq/chizer-1.0) < CHI_SQR_LIMIT) and  (test < TEST_LIMIT):
618            print ' Convergence achieved.'
619            return chisq,f,image_to_data(f, G)     # solution FOUND returns here
620    print ' No convergence! Try increasing Error multiplier.'
621    return chisq,f,image_to_data(f, G)       # no solution after IterMax iterations
622
623   
624################################################################################
625#### MaxEnt testing stuff
626################################################################################
627
628def print_vec(text, a):
629    '''print the contents of a vector to the console'''
630    n = a.shape[0]
631    print "%s[ = (" % text,
632    for i in range(n):
633        s = " %g, " % a[i]
634        print s,
635    print ")"
636
637def print_arr(text, a):
638    '''print the contents of an array to the console'''
639    n, m = a.shape
640    print "%s[][] = (" % text
641    for i in range(n):
642        print " (",
643        for j in range(m):
644            print " %g, " % a[i][j],
645        print "),"
646    print ")"
647
648def test_MaxEnt_SB(report=True):
649    def readTextData(filename):
650        '''return q, I, dI from a 3-column text file'''
651        if not os.path.exists(filename):
652            raise Exception("file not found: " + filename)
653        buf = [line.split() for line in open(filename, 'r').readlines()]
654        M = len(buf)
655        buf = zip(*buf)         # transpose rows and columns
656        q  = numpy.array(buf[0], dtype=numpy.float64)
657        I  = numpy.array(buf[1], dtype=numpy.float64)
658        dI = numpy.array(buf[2], dtype=numpy.float64)
659        return q, I, dI
660    print "MaxEnt_SB: "
661    test_data_file = os.path.join( 'testinp', 'test.sas')
662    rhosq = 100     # scattering contrast, 10^20 1/cm^-4
663    bkg   = 0.1     #   I = I - bkg
664    dMin, dMax, nRadii = 25, 9000, 40
665    defaultDistLevel = 1.0e-6
666    IterMax = 40
667    errFac = 1.05
668   
669    r    = numpy.logspace(math.log10(dMin), math.log10(dMax), nRadii)/2
670    dr   = r * (r[1]/r[0] - 1)          # step size
671    f_dr = numpy.ndarray((nRadii)) * 0  # volume fraction histogram
672    b    = numpy.ndarray((nRadii)) * 0 + defaultDistLevel  # MaxEnt "sky background"
673   
674    qVec, I, dI = readTextData(test_data_file)
675    G = G_matrix(qVec,r,rhosq,SphereFF,SphereVol,args=())
676   
677    chisq,f_dr,Ic = MaxEnt_SB(I - bkg, dI*errFac, b, IterMax, G, report=report)
678    if f_dr is None:
679        print "no solution"
680        return
681   
682    print "solution reached"
683    for a,b,c in zip(r.tolist(), dr.tolist(), f_dr.tolist()):
684        print '%10.4f %10.4f %12.4g'%(a,b,c)
685
686def tests():
687    test_MaxEnt_SB(report=True)
688
689if __name__ == '__main__':
690    tests()
691   
692###############################################################################
693#### Size distribution
694###############################################################################
695
696def SizeDistribution(Profile,ProfDict,Limits,Substances,Sample,data):
697    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
698        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
699        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
700        'Unified disk':[UniDiskFF,UniDiskVol]}
701    Shape = data['Size']['Shape'][0]
702    Parms = data['Size']['Shape'][1:]
703    if data['Size']['logBins']:
704        Bins = np.logspace(np.log10(data['Size']['MinDiam']),np.log10(data['Size']['MaxDiam']),
705            data['Size']['Nbins']+1,True)/2.        #make radii
706    else:
707        Bins = np.linspace(data['Size']['MinDiam'],data['Size']['MaxDiam'],
708            data['Size']['Nbins']+1,True)/2.        #make radii
709    Dbins = np.diff(Bins)
710    Bins = Bins[:-1]+Dbins/2.
711    BinsBack = np.ones_like(Bins)*1.e-6
712    Back = data['Back']
713    Q,Io,wt,Ic,Ib = Profile[:5]
714    Qmin = Limits[1][0]
715    Qmax = Limits[1][1]
716    Contrast = Sample['Contrast'][1]
717    wtFactor = ProfDict['wtFactor']
718    Ibeg = np.searchsorted(Q,Qmin)
719    Ifin = np.searchsorted(Q,Qmax)
720    if Back[1]:
721        Ib = Back[0]
722        Ic[Ibeg:Ifin] = Back[0]
723    Gmat = G_matrix(Q[Ibeg:Ifin],Bins,Contrast,shapes[Shape][0],shapes[Shape][1],args=Parms)
724    chisq,BinMag,Ic[Ibeg:Ifin] = MaxEnt_SB(Io[Ibeg:Ifin]-Back[0],1./np.sqrt(wtFactor*wt[Ibeg:Ifin]),BinsBack,
725        data['Size']['MaxEnt']['Niter'],Gmat,report=True)
726    print ' Final chi^2: %.3f'%(chisq)
727    Vols = shapes[Shape][1](Bins,Parms)
728    data['Size']['Distribution'] = [Bins,Dbins,BinMag/(2.*Dbins)]
729       
730   
Note: See TracBrowser for help on using the repository browser.