source: trunk/GSASIIsasd.py @ 1249

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

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