source: trunk/GSASIIsasd.py @ 1244

Last change on this file since 1244 was 1244, checked in by vondreele, 9 years ago

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