source: trunk/GSASIIsasd.py @ 1302

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

merge messup fixed!
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 49.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-04-25 17:56:35 +0000 (Fri, 25 Apr 2014) $
10# $Author: vondreele $
11# $Revision: 1302 $
12# $URL: trunk/GSASIIsasd.py $
13# $Id: GSASIIsasd.py 1302 2014-04-25 17:56:35Z 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#import pdb
29
30import GSASIIpath
31GSASIIpath.SetVersionNumber("$Revision: 1302 $")
32import GSASIIlattice as G2lat
33import GSASIIspc as G2spc
34import GSASIIElem as G2elem
35import GSASIIgrid as G2gd
36import GSASIIIO as G2IO
37import GSASIImath as G2mth
38import GSASIIpwd as G2pwd
39
40# trig functions in degrees
41sind = lambda x: math.sin(x*math.pi/180.)
42asind = lambda x: 180.*math.asin(x)/math.pi
43tand = lambda x: math.tan(x*math.pi/180.)
44atand = lambda x: 180.*math.atan(x)/math.pi
45atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
46cosd = lambda x: math.cos(x*math.pi/180.)
47acosd = lambda x: 180.*math.acos(x)/math.pi
48rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
49#numpy versions
50npsind = lambda x: np.sin(x*np.pi/180.)
51npasind = lambda x: 180.*np.arcsin(x)/math.pi
52npcosd = lambda x: np.cos(x*math.pi/180.)
53npacosd = lambda x: 180.*np.arccos(x)/math.pi
54nptand = lambda x: np.tan(x*math.pi/180.)
55npatand = lambda x: 180.*np.arctan(x)/np.pi
56npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
57npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave
58npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
59   
60###############################################################################
61#### Particle form factors
62###############################################################################
63
64def SphereFF(Q,R,args=()):
65    ''' Compute hard sphere form factor - can use numpy arrays
66    param float Q: Q value array (usually in A-1)
67    param float R: sphere radius (Usually in A - must match Q-1 units)
68    param array args: ignored
69    returns float: form factors as array as needed
70    '''
71    QR = Q[:,np.newaxis]*R
72    return (3./(QR**3))*(np.sin(QR)-(QR*np.cos(QR)))
73   
74def SpheroidFF(Q,R,args):
75    ''' Compute form factor of cylindrically symmetric ellipsoid (spheroid)
76    - can use numpy arrays for R & AR; will return corresponding numpy array
77    param float Q : Q value array (usually in A-1)
78    param float R: radius along 2 axes of spheroid
79    param array args: [float AR]: aspect ratio so 3rd axis = R*AR
80    returns float: form factors as array as needed
81    '''
82    NP = 50
83    AR = args[0]
84    if 0.99 < AR < 1.01:
85        return SphereFF(Q,R,0)
86    else:
87        cth = np.linspace(0,1.,NP)
88        Rct = R[:,np.newaxis]*np.sqrt(1.+(AR**2-1.)*cth**2)
89        return np.sqrt(np.sum(SphereFF(Q[:,np.newaxis],Rct,0)**2,axis=2)/NP)
90           
91def CylinderFF(Q,R,args):
92    ''' Compute form factor for cylinders - can use numpy arrays
93    param float Q: Q value array (A-1)
94    param float R: cylinder radius (A)
95    param array args: [float L]: cylinder length (A)
96    returns float: form factor
97    '''
98    L = args[0]
99    NP = 200
100    alp = np.linspace(0,np.pi/2.,NP)
101    QL = Q[:,np.newaxis]*L
102    LBessArg = 0.5*QL[:,:,np.newaxis]**np.cos(alp)
103    LBess = np.where(LBessArg<1.e-6,1.,np.sin(LBessArg)/LBessArg)
104    QR = Q[:,np.newaxis]*R
105    SBessArg = QR[:,:,np.newaxis]*np.sin(alp)
106    SBess = np.where(SBessArg<1.e-6,0.5,scsp.jv(1,SBessArg)/SBessArg)
107    LSBess = LBess*SBess
108    return np.sqrt(2.*np.pi*np.sum(np.sin(alp)*LSBess**2,axis=2)/NP)
109   
110def CylinderDFF(Q,L,args):
111    ''' Compute form factor for cylinders - can use numpy arrays
112    param float Q: Q value array (A-1)
113    param float L: cylinder half length (A)
114    param array args: [float R]: cylinder radius (A)
115    returns float: form factor
116    '''
117    R = args[0]
118    return CylinderFF(Q,R,args=[2.*L])   
119   
120def CylinderARFF(Q,R,args): 
121    ''' Compute form factor for cylinders - can use numpy arrays
122    param float Q: Q value array (A-1)
123    param float R: cylinder radius (A)
124    param array args: [float AR]: cylinder aspect ratio = L/D = L/2R
125    returns float: form factor
126    '''
127    AR = args[0]
128    return CylinderFF(Q,R,args=[2.*R*AR])   
129   
130def UniSphereFF(Q,R,args=0):
131    ''' Compute form factor for unified sphere - can use numpy arrays
132    param float Q: Q value array (A-1)
133    param float R: cylinder radius (A)
134    param array args: ignored
135    returns float: form factor
136    '''
137    Rg = np.sqrt(3./5.)*R
138    B = np.pi*1.62/(Rg**4)    #are we missing *np.pi? 1.62 = 6*(3/5)**2/(4/3) sense?
139    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg/np.sqrt(6)))**3
140    return np.sqrt(np.exp((-Q[:,np.newaxis]**2*Rg**2)/3.)+(B/QstV**4))
141   
142def UniRodFF(Q,R,args):
143    ''' Compute form factor for unified rod - can use numpy arrays
144    param float Q: Q value array (A-1)
145    param float R: cylinder radius (A)
146    param array args: [float R]: cylinder radius (A)
147    returns float: form factor
148    '''
149    L = args[0]
150    Rg2 = np.sqrt(R**2/2+L**2/12)
151    B2 = np.pi/L
152    Rg1 = np.sqrt(3.)*R/2.
153    G1 = (2./3.)*R/L
154    B1 = 4.*(L+R)/(R**3*L**2)
155    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg2/np.sqrt(6)))**3
156    FF = np.exp(-Q[:,np.newaxis]**2*Rg2**2/3.)+(B2/QstV)*np.exp(-Rg1**2*Q[:,np.newaxis]**2/3.)
157    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg1/np.sqrt(6)))**3
158    FF += G1*np.exp(-Q[:,np.newaxis]**2*Rg1**2/3.)+(B1/QstV**4)
159    return np.sqrt(FF)
160   
161def UniRodARFF(Q,R,args):
162    ''' Compute form factor for unified rod of fixed aspect ratio - can use numpy arrays
163    param float Q: Q value array (A-1)
164    param float R: cylinder radius (A)
165    param array args: [float AR]: cylinder aspect ratio = L/D = L/2R
166    returns float: form factor
167    '''
168    AR = args[0]
169    return UniRodFF(Q,R,args=[2.*AR*R,])
170   
171def UniDiskFF(Q,R,args):
172    ''' Compute form factor for unified disk - can use numpy arrays
173    param float Q: Q value array (A-1)
174    param float R: cylinder radius (A)
175    param array args: [float T]: disk thickness (A)
176    returns float: form factor
177    '''
178    T = args[0]
179    Rg2 = np.sqrt(R**2/2.+T**2/12.)
180    B2 = 2./R**2
181    Rg1 = np.sqrt(3.)*T/2.
182    RgC2 = 1.1*Rg1
183    G1 = (2./3.)*(T/R)**2
184    B1 = 4.*(T+R)/(R**3*T**2)
185    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg2/np.sqrt(6)))**3
186    FF = np.exp(-Q[:,np.newaxis]**2*Rg2**2/3.)+(B2/QstV**2)*np.exp(-RgC2**2*Q[:,np.newaxis]**2/3.)
187    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg1/np.sqrt(6)))**3
188    FF += G1*np.exp(-Q[:,np.newaxis]**2*Rg1**2/3.)+(B1/QstV**4)
189    return np.sqrt(FF)
190   
191def UniTubeFF(Q,R,args):
192    ''' Compute form factor for unified tube - can use numpy arrays
193    assumes that core of tube is same as the matrix/solvent so contrast
194    is from tube wall vs matrix
195    param float Q: Q value array (A-1)
196    param float R: cylinder radius (A)
197    param array args: [float L,T]: tube length & wall thickness(A)
198    returns float: form factor
199    '''
200    L,T = args[:2]
201    Ri = R-T
202    DR2 = R**2-Ri**2
203    Vt = np.pi*DR2*L
204    Rg3 = np.sqrt(DR2/2.+L**2/12.)
205    B1 = 4.*np.pi**2*(DR2+L*(R+Ri))/Vt**2
206    B2 = np.pi**2*T/Vt
207    B3 = np.pi/L
208    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg3/np.sqrt(6)))**3
209    FF = np.exp(-Q[:,np.newaxis]**2*Rg3**2/3.)+(B3/QstV)*np.exp(-Q[:,np.newaxis]**2*R**2/3.)
210    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*R/np.sqrt(6)))**3
211    FF += (B2/QstV**2)*np.exp(-Q[:,np.newaxis]**2*T**2/3.)
212    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*T/np.sqrt(6)))**3
213    FF += B1/QstV**4
214    return np.sqrt(FF)
215
216###############################################################################
217#### Particle volumes
218###############################################################################
219
220def SphereVol(R,args=()):
221    ''' Compute volume of sphere
222    - numpy array friendly
223    param float R: sphere radius
224    param array args: ignored
225    returns float: volume
226    '''
227    return (4./3.)*np.pi*R**3
228
229def SpheroidVol(R,args):
230    ''' Compute volume of cylindrically symmetric ellipsoid (spheroid)
231    - numpy array friendly
232    param float R: radius along 2 axes of spheroid
233    param array args: [float AR]: aspect ratio so radius of 3rd axis = R*AR
234    returns float: volume
235    '''
236    AR = args[0]
237    return AR*SphereVol(R)
238   
239def CylinderVol(R,args):
240    ''' Compute cylinder volume for radius & length
241    - numpy array friendly
242    param float R: diameter (A)
243    param array args: [float L]: length (A)
244    returns float:volume (A^3)
245    '''
246    L = args[0]
247    return np.pi*L*R**2
248   
249def CylinderDVol(L,args):
250    ''' Compute cylinder volume for length & diameter
251    - numpy array friendly
252    param float: L half length (A)
253    param array args: [float D]: diameter (A)
254    returns float:volume (A^3)
255    '''
256    D = args[0]
257    return CylinderVol(D/2.,[2.*L,])
258   
259def CylinderARVol(R,args):
260    ''' Compute cylinder volume for radius & aspect ratio = L/D
261    - numpy array friendly
262    param float: R radius (A)
263    param array args: [float AR]: =L/D=L/2R aspect ratio
264    returns float:volume
265    '''
266    AR = args[0]
267    return CylinderVol(R,[2.*R*AR,])
268   
269def UniSphereVol(R,args=()):
270    ''' Compute volume of sphere
271    - numpy array friendly
272    param float R: sphere radius
273    param array args: ignored
274    returns float: volume
275    '''
276    return SphereVol(R)
277   
278def UniRodVol(R,args):
279    ''' Compute cylinder volume for radius & length
280    - numpy array friendly
281    param float R: diameter (A)
282    param array args: [float L]: length (A)
283    returns float:volume (A^3)
284    '''
285    L = args[0]
286    return CylinderVol(R,[L,])
287   
288def UniRodARVol(R,args):
289    ''' Compute rod volume for radius & aspect ratio
290    - numpy array friendly
291    param float R: diameter (A)
292    param array args: [float AR]: =L/D=L/2R aspect ratio
293    returns float:volume (A^3)
294    '''
295    AR = args[0]
296    return CylinderARVol(R,[AR,])
297   
298def UniDiskVol(R,args):
299    ''' Compute disk volume for radius & thickness
300    - numpy array friendly
301    param float R: diameter (A)
302    param array args: [float T]: thickness
303    returns float:volume (A^3)
304    '''
305    T = args[0]
306    return CylinderVol(R,[T,])
307   
308def UniTubeVol(R,args):
309    ''' Compute tube volume for radius, length & wall thickness
310    - numpy array friendly
311    param float R: diameter (A)
312    param array args: [float L,T]: tube length & wall thickness(A)
313    returns float: volume (A^3) of tube wall
314    '''
315    L,T = args[:2]
316    return CylinderVol(R,[L,])-CylinderVol(R-T,[L,])
317   
318################################################################################
319#### Distribution functions & their cumulative fxns
320################################################################################
321
322def LogNormalDist(x,pos,args):
323    ''' Standard LogNormal distribution - numpy friendly on x axis
324    ref: http://www.itl.nist.gov/div898/handbook/index.htm 1.3.6.6.9
325    param float x: independent axis (can be numpy array)
326    param float pos: location of distribution
327    param float scale: width of distribution (m)
328    param float shape: shape - (sigma of log(LogNormal))
329    returns float: LogNormal distribution
330    '''
331    scale,shape = args
332    return np.exp(-np.log((x-pos)/scale)**2/(2.*shape**2))/(np.sqrt(2.*np.pi)*(x-pos)*shape)
333   
334def GaussDist(x,pos,args):
335    ''' Standard Normal distribution - numpy friendly on x axis
336    param float x: independent axis (can be numpy array)
337    param float pos: location of distribution
338    param float scale: width of distribution (sigma)
339    param float shape: not used
340    returns float: Normal distribution
341    '''
342    scale = args[0]
343    return (1./(scale*np.sqrt(2.*np.pi)))*np.exp(-(x-pos)**2/(2.*scale**2))
344   
345def LSWDist(x,pos,args=[]):
346    ''' Lifshitz-Slyozov-Wagner Ostwald ripening distribution - numpy friendly on x axis
347    ref:
348    param float x: independent axis (can be numpy array)
349    param float pos: location of distribution
350    param float scale: not used
351    param float shape: not used
352    returns float: LSW distribution   
353    '''
354    redX = x/pos
355    result = (81.*2**(-5/3.))*(redX**2*np.exp(-redX/(1.5-redX)))/((1.5-redX)**(11/3.)*(3.-redX)**(7/3.))
356   
357    return np.nan_to_num(result/pos)
358   
359def SchulzZimmDist(x,pos,args):
360    ''' Schulz-Zimm macromolecule distribution - numpy friendly on x axis
361    ref: http://goldbook.iupac.org/S05502.html
362    param float x: independent axis (can be numpy array)
363    param float pos: location of distribution
364    param float scale: width of distribution (sigma)
365    param float shape: not used
366    returns float: Schulz-Zimm distribution
367    '''
368    scale = args[0]
369    b = (2.*pos/scale)**2
370    a = b/pos
371    if b<70:    #why bother?
372        return (a**(b+1.))*x**b*np.exp(-a*x)/scsp.gamma(b+1.)
373    else:
374        return np.exp((b+1.)*np.log(a)-scsp.gammaln(b+1.)+b*np.log(x)-(a*x))
375           
376def LogNormalCume(x,pos,args):
377    ''' Standard LogNormal cumulative distribution - numpy friendly on x axis
378    ref: http://www.itl.nist.gov/div898/handbook/index.htm 1.3.6.6.9
379    param float x: independent axis (can be numpy array)
380    param float pos: location of distribution
381    param float scale: width of distribution (sigma)
382    param float shape: shape parameter
383    returns float: LogNormal cumulative distribution
384    '''
385    scale,shape = args
386    lredX = np.log((x-pos)/scale)
387    return (scsp.erf((lredX/shape)/np.sqrt(2.))+1.)/2.
388   
389def GaussCume(x,pos,args):
390    ''' Standard Normal cumulative distribution - numpy friendly on x axis
391    param float x: independent axis (can be numpy array)
392    param float pos: location of distribution
393    param float scale: width of distribution (sigma)
394    param float shape: not used
395    returns float: Normal cumulative distribution
396    '''
397    scale = args[0]
398    redX = (x-pos)/scale
399    return (scsp.erf(redX/np.sqrt(2.))+1.)/2.
400   
401def LSWCume(x,pos,args=[]):
402    ''' Lifshitz-Slyozov-Wagner Ostwald ripening cumulative distribution - numpy friendly on x axis
403    param float x: independent axis (can be numpy array)
404    param float pos: location of distribution
405    param float scale: not used
406    param float shape: not used
407    returns float: LSW cumulative distribution
408    '''
409    nP = 500
410    xMin,xMax = [0.,2.*pos]
411    X = np.linspace(xMin,xMax,nP,True)
412    fxn = LSWDist(X,pos)
413    mat = np.outer(np.ones(nP),fxn)
414    cume = np.sum(np.tril(mat),axis=1)/np.sum(fxn)
415    return np.interp(x,X,cume,0,1)
416   
417def SchulzZimmCume(x,pos,args):
418    ''' Schulz-Zimm cumulative distribution - numpy friendly on x axis
419    param float x: independent axis (can be numpy array)
420    param float pos: location of distribution
421    param float scale: width of distribution (sigma)
422    param float shape: not used
423    returns float: Normal distribution
424    '''
425    scale = args[0]
426    nP = 500
427    xMin = np.max([0.,pos-4.*scale])
428    xMax = np.min([pos+4.*scale,1.e5])
429    X = np.linspace(xMin,xMax,nP,True)
430    fxn = LSWDist(X,pos)
431    mat = np.outer(np.ones(nP),fxn)
432    cume = np.sum(np.tril(mat),axis=1)/np.sum(fxn)
433    return np.interp(x,X,cume,0,1)
434   
435    return []
436   
437         
438################################################################################
439##### SB-MaxEnt
440################################################################################
441
442def G_matrix(q,r,contrast,FFfxn,Volfxn,args=()):
443    '''Calculates the response matrix :math:`G(Q,r)`
444   
445    :param float q: :math:`Q`
446    :param float r: :math:`r`
447    :param float contrast: :math:`|\\Delta\\rho|^2`, the scattering contrast
448    :param function FFfxn: form factor function FF(q,r,args)
449    :param function Volfxn: volume function Vol(r,args)
450    :returns float: G(Q,r)
451    '''
452    FF = FFfxn(q,r,args)
453    Vol = Volfxn(r,args)
454    return 1.e-4*(contrast*Vol*FF**2).T     #10^-20 vs 10^-24
455   
456'''
457sbmaxent
458
459Entropy maximization routine as described in the article
460J Skilling and RK Bryan; MNRAS 211 (1984) 111 - 124.
461("MNRAS": "Monthly Notices of the Royal Astronomical Society")
462
463:license: Copyright (c) 2013, UChicago Argonne, LLC
464:license: This file is distributed subject to a Software License Agreement found
465     in the file LICENSE that is included with this distribution.
466
467References:
468
4691. J Skilling and RK Bryan; MON NOT R ASTR SOC 211 (1984) 111 - 124.
4702. JA Potton, GJ Daniell, and BD Rainford; Proc. Workshop
471   Neutron Scattering Data Analysis, Rutherford
472   Appleton Laboratory, UK, 1986; ed. MW Johnson,
473   IOP Conference Series 81 (1986) 81 - 86, Institute
474   of Physics, Bristol, UK.
4753. ID Culverwell and GP Clarke; Ibid. 87 - 96.
4764. JA Potton, GK Daniell, & BD Rainford,
477   J APPL CRYST 21 (1988) 663 - 668.
4785. JA Potton, GJ Daniell, & BD Rainford,
479   J APPL CRYST 21 (1988) 891 - 897.
480
481'''
482
483class MaxEntException(Exception): 
484    '''Any exception from this module'''
485    pass
486
487def MaxEnt_SB(datum, sigma, G, base, IterMax, image_to_data=None, data_to_image=None, report=False):
488    '''
489    do the complete Maximum Entropy algorithm of Skilling and Bryan
490   
491    :param float datum[]:
492    :param float sigma[]:
493    :param float[][] G: transformation matrix
494    :param float base[]:
495    :param int IterMax:
496    :param obj image_to_data: opus function (defaults to opus)
497    :param obj data_to_image: tropus function (defaults to tropus)
498   
499    :returns float[]: :math:`f(r) dr`
500    '''
501       
502    TEST_LIMIT        = 0.05                    # for convergence
503    CHI_SQR_LIMIT     = 0.01                    # maximum difference in ChiSqr for a solution
504    SEARCH_DIRECTIONS = 3                       # <10.  This code requires value = 3
505    RESET_STRAYS      = 1                       # was 0.001, correction of stray negative values
506    DISTANCE_LIMIT_FACTOR = 0.1                 # limitation on df to constrain runaways
507   
508    MAX_MOVE_LOOPS    = 500                     # for no solution in routine: move,
509    MOVE_PASSES       = 0.001                   # convergence test in routine: move
510
511    def tropus (data, G):
512        '''
513        tropus: transform data-space -> solution-space:  [G] * data
514       
515        default definition, caller can use this definition or provide an alternative
516       
517        :param float[M] data: observations, ndarray of shape (M)
518        :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
519        :returns float[N]: calculated image, ndarray of shape (N)
520        '''
521        return G.dot(data)
522
523    def opus (image, G):
524        '''
525        opus: transform solution-space -> data-space:  [G]^tr * image
526       
527        default definition, caller can use this definition or provide an alternative
528       
529        :param float[N] image: solution, ndarray of shape (N)
530        :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
531        :returns float[M]: calculated data, ndarray of shape (M)
532        '''
533        return np.dot(G.T,image)    #G.transpose().dot(image)
534
535    def Dist(s2, beta):
536        '''measure the distance of this possible solution'''
537        w = 0
538        n = beta.shape[0]
539        for k in range(n):
540            z = -sum(s2[k] * beta)
541            w += beta[k] * z
542        return w
543   
544    def ChiNow(ax, c1, c2, s1, s2):
545        '''
546        ChiNow
547       
548        :returns tuple: (ChiNow computation of ``w``, beta)
549        '''
550       
551        bx = 1 - ax
552        a =   bx * c2  -  ax * s2
553        b = -(bx * c1  -  ax * s1)
554   
555        beta = ChoSol(a, b)
556        w = 1.0
557        for k in range(SEARCH_DIRECTIONS):
558            w += beta[k] * (c1[k] + 0.5*sum(c2[k] * beta))
559        return w, beta
560   
561    def ChoSol(a, b):
562        '''
563        ChoSol: ? chop the solution vectors ?
564       
565        :returns: new vector beta
566        '''
567        n = b.shape[0]
568        fl = np.zeros((n,n))
569        bl = np.zeros_like(b)
570       
571        #print_arr("ChoSol: a", a)
572        #print_vec("ChoSol: b", b)
573   
574        if (a[0][0] <= 0):
575            msg = "ChoSol: a[0][0] = " 
576            msg += str(a[0][0])
577            msg += '  Value must be positive'
578            raise MaxEntException(msg)
579   
580        # first, compute fl from a
581        # note fl is a lower triangular matrix
582        fl[0][0] = math.sqrt (a[0][0])
583        for i in (1, 2):
584            fl[i][0] = a[i][0] / fl[0][0]
585            for j in range(1, i+1):
586                z = 0.0
587                for k in range(j):
588                    z += fl[i][k] * fl[j][k]
589                    #print "ChoSol: %d %d %d  z = %lg" % ( i, j, k, z)
590                z = a[i][j] - z
591                if j == i:
592                    y = math.sqrt(z)
593                else:
594                    y = z / fl[j][j]
595                fl[i][j] = y
596        #print_arr("ChoSol: fl", fl)
597   
598        # next, compute bl from fl and b
599        bl[0] = b[0] / fl[0][0]
600        for i in (1, 2):
601            z = 0.0
602            for k in range(i):
603                z += fl[i][k] * bl[k]
604                #print "\t", i, k, z
605            bl[i] = (b[i] - z) / fl[i][i]
606        #print_vec("ChoSol: bl", bl)
607   
608        # last, compute beta from bl and fl
609        beta = np.empty((n))
610        beta[-1] = bl[-1] / fl[-1][-1]
611        for i in (1, 0):
612            z = 0.0
613            for k in range(i+1, n):
614                z += fl[k][i] * beta[k]
615                #print "\t\t", i, k, 'z=', z
616            beta[i] = (bl[i] - z) / fl[i][i]
617        #print_vec("ChoSol: beta", beta)
618   
619        return beta
620
621    def MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2):
622        '''
623        move beta one step closer towards the solution
624        '''
625        a_lower, a_upper = 0., 1.          # bracket  "a"
626        cmin, beta = ChiNow (a_lower, c1, c2, s1, s2)
627        #print "MaxEntMove: cmin = %g" % cmin
628        if cmin*chisq > chizer:
629            ctarg = (1.0 + cmin)/2
630        else:
631            ctarg = chizer/chisq
632        f_lower = cmin - ctarg
633        c_upper, beta = ChiNow (a_upper, c1, c2, s1, s2)
634        f_upper = c_upper - ctarg
635   
636        fx = 2*MOVE_PASSES      # just to start off
637        loop = 1
638        while abs(fx) >= MOVE_PASSES and loop <= MAX_MOVE_LOOPS:
639            a_new = (a_lower + a_upper) * 0.5           # search by bisection
640            c_new, beta = ChiNow (a_new, c1, c2, s1, s2)
641            fx = c_new - ctarg
642            # tighten the search range for the next pass
643            if f_lower*fx > 0:
644                a_lower, f_lower = a_new, fx
645            if f_upper*fx > 0:
646                a_upper, f_upper = a_new, fx
647            loop += 1
648   
649        if abs(fx) >= MOVE_PASSES or loop > MAX_MOVE_LOOPS:
650            msg = "MaxEntMove: Loop counter = " 
651            msg += str(MAX_MOVE_LOOPS)
652            msg += '  No convergence in alpha chop'
653            raise MaxEntException(msg)
654   
655        w = Dist (s2, beta);
656        m = SEARCH_DIRECTIONS
657        if (w > DISTANCE_LIMIT_FACTOR*fSum/blank):        # invoke the distance penalty, SB eq. 17
658            for k in range(m):
659                beta[k] *= math.sqrt (fSum/(blank*w))
660        chtarg = ctarg * chisq
661        return w, chtarg, loop, a_new, fx, beta
662       
663#MaxEnt_SB starts here   
664
665    if image_to_data == None:
666        image_to_data = opus
667    if data_to_image == None:
668        data_to_image = tropus
669    n   = len(base)
670    npt = len(datum)
671
672    # Note that the order of subscripts for
673    # "xi" and "eta" has been reversed from
674    # the convention used in the FORTRAN version
675    # to enable parts of them to be passed as
676    # as vectors to "image_to_data" and "data_to_image".
677    xi      = np.zeros((SEARCH_DIRECTIONS, n))
678    eta     = np.zeros((SEARCH_DIRECTIONS, npt))
679    beta    = np.zeros((SEARCH_DIRECTIONS))
680    # s1      = np.zeros((SEARCH_DIRECTIONS))
681    # c1      = np.zeros((SEARCH_DIRECTIONS))
682    s2      = np.zeros((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
683    c2      = np.zeros((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
684
685    # TODO: replace blank (scalar) with base (vector)
686    blank = sum(base) / len(base)   # use the average value of base
687
688    chizer, chtarg = npt*1.0, npt*1.0
689    f = base * 1.0                              # starting distribution is base
690
691    fSum  = sum(f)                              # find the sum of the f-vector
692    z = (datum - image_to_data (f, G)) / sigma  # standardized residuals, SB eq. 3
693    chisq = sum(z*z)                            # Chi^2, SB eq. 4
694
695    for iter in range(IterMax):
696        ox = -2 * z / sigma                        # gradient of Chi^2
697
698        cgrad = data_to_image (ox, G)              # cgrad[i] = del(C)/del(f[i]), SB eq. 8
699        sgrad = -np.log(f/base) / (blank*math.exp (1.0))  # sgrad[i] = del(S)/del(f[i])
700        snorm = math.sqrt(sum(f * sgrad*sgrad))    # entropy term, SB eq. 22
701        cnorm = math.sqrt(sum(f * cgrad*cgrad))    # ChiSqr term, SB eq. 22
702        tnorm = sum(f * sgrad * cgrad)             # norm for gradient term TEST
703
704        a = 1.0
705        b = 1.0 / cnorm
706        if iter == 0:
707            test = 0.0     # mismatch between entropy and ChiSquared gradients
708        else:
709            test = math.sqrt( ( 1.0 - tnorm/(snorm*cnorm) )/2 ) # SB eq. 37?
710            a = 0.5 / (snorm * test)
711            b *= 0.5 / test
712        xi[0] = f * cgrad / cnorm
713        xi[1] = f * (a * sgrad - b * cgrad)
714
715        eta[0] = image_to_data (xi[0], G);          # image --> data
716        eta[1] = image_to_data (xi[1], G);          # image --> data
717        ox = eta[1] / (sigma * sigma)
718        xi[2] = data_to_image (ox, G);              # data --> image
719        a = 1.0 / math.sqrt(sum(f * xi[2]*xi[2]))
720        xi[2] = f * xi[2] * a
721        eta[2] = image_to_data (xi[2], G)           # image --> data
722       
723#         print_arr("MaxEnt: eta.transpose()", eta.transpose())
724#         print_arr("MaxEnt: xi.transpose()", xi.transpose())
725
726        # prepare the search directions for the conjugate gradient technique
727        c1 = xi.dot(cgrad) / chisq                          # C_mu, SB eq. 24
728        s1 = xi.dot(sgrad)                          # S_mu, SB eq. 24
729#         print_vec("MaxEnt: c1", c1)
730#         print_vec("MaxEnt: s1", s1)
731
732        for k in range(SEARCH_DIRECTIONS):
733            for l in range(k+1):
734                c2[k][l] = 2 * sum(eta[k] * eta[l] / sigma/sigma) / chisq
735                s2[k][l] = -sum(xi[k] * xi[l] / f) / blank
736#         print_arr("MaxEnt: c2", c2)
737#         print_arr("MaxEnt: s2", s2)
738
739        # reflect across the body diagonal
740        for k, l in ((0,1), (0,2), (1,2)):
741            c2[k][l] = c2[l][k]                     #  M_(mu,nu)
742            s2[k][l] = s2[l][k]                     #  g_(mu,nu)
743 
744        beta[0] = -0.5 * c1[0] / c2[0][0]
745        beta[1] = 0.0
746        beta[2] = 0.0
747        if (iter > 0):
748            w, chtarg, loop, a_new, fx, beta = MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2)
749
750        f_old = f.copy()    # preserve the last image
751        f += xi.transpose().dot(beta)   # move the image towards the solution, SB eq. 25
752       
753        # As mentioned at the top of p.119,
754        # need to protect against stray negative values.
755        # In this case, set them to RESET_STRAYS * base[i]
756        #f = f.clip(RESET_STRAYS * blank, f.max())
757        for i in range(n):
758            if f[i] <= 0.0:
759                f[i] = RESET_STRAYS * base[i]
760        df = f - f_old
761        fSum = sum(f)
762        fChange = sum(df)
763
764        # calculate the normalized entropy
765        S = sum((f/fSum) * np.log(f/fSum))      # normalized entropy, S&B eq. 1
766        z = (datum - image_to_data (f, G)) / sigma  # standardized residuals
767        chisq = sum(z*z)                            # report this ChiSq
768
769        if report:
770            print " MaxEnt trial/max: %3d/%3d" % ((iter+1), IterMax)
771            print " Residual: %5.2lf%% Entropy: %8lg" % (100*test, S)
772            if iter > 0:
773                value = 100*( math.sqrt(chisq/chtarg)-1)
774            else:
775                value = 0
776      #      print " %12.5lg %10.4lf" % ( math.sqrt(chtarg/npt), value )
777            print " Function sum: %.6lg Change from last: %.2lf%%\n" % (fSum,100*fChange/fSum)
778
779        # See if we have finished our task.
780        # do the hardest test first
781        if (abs(chisq/chizer-1.0) < CHI_SQR_LIMIT) and  (test < TEST_LIMIT):
782            print ' Convergence achieved.'
783            return chisq,f,image_to_data(f, G)     # solution FOUND returns here
784    print ' No convergence! Try increasing Error multiplier.'
785    return chisq,f,image_to_data(f, G)       # no solution after IterMax iterations
786
787   
788###############################################################################
789#### IPG/TNNLS Routines
790###############################################################################
791
792def IPG(datum,sigma,G,Bins,Dbins,IterMax,Qvec=[],approach=0.8,Power=-1,report=False):
793    ''' An implementation of the Interior-Point Gradient method of
794    Michael Merritt & Yin Zhang, Technical Report TR04-08, Dept. of Comp. and
795    Appl. Math., Rice Univ., Houston, Texas 77005, U.S.A. found on the web at
796    http://www.caam.rice.edu/caam/trs/2004/TR04-08.pdf
797    Problem addressed: Total Non-Negative Least Squares (TNNLS)
798    :param float datum[]:
799    :param float sigma[]:
800    :param float[][] G: transformation matrix
801    :param int IterMax:
802    :param float Qvec: data positions for Power = 0-4
803    :param float approach: 0.8 default fitting parameter
804    :param int Power: 0-4 for Q^Power weighting, -1 to use input sigma
805   
806    '''
807    if Power < 0: 
808        GmatE = G/sigma[:np.newaxis]
809        IntE = datum/sigma
810        pwr = 0
811        QvecP = np.ones_like(datum)
812    else:
813        GmatE = G[:]
814        IntE = datum[:]
815        pwr = Power
816        QvecP = Qvec**pwr
817    Amat = GmatE*QvecP[:np.newaxis]
818    AAmat = np.inner(Amat,Amat)
819    Bvec = IntE*QvecP
820    Xw = np.ones_like(Bins)*1.e-32
821    calc = np.dot(G.T,Xw)
822    nIter = 0
823    err = 10.
824    while (nIter<IterMax) and (err > 1.):
825        #Step 1 in M&Z paper:
826        Qk = np.inner(AAmat,Xw)-np.inner(Amat,Bvec)
827        Dk = Xw/np.inner(AAmat,Xw)
828        Pk = -Dk*Qk
829        #Step 2 in M&Z paper:
830        alpSt = -np.inner(Pk,Qk)/np.inner(Pk,np.inner(AAmat,Pk))
831        alpWv = -Xw/Pk
832        AkSt = [np.where(Pk[i]<0,np.min((approach*alpWv[i],alpSt)),alpSt) for i in range(Pk.shape[0])]
833        #Step 3 in M&Z paper:
834        shift = AkSt*Pk
835        Xw += shift
836        #done IPG; now results
837        nIter += 1
838        calc = np.dot(G.T,Xw)
839        chisq = np.sum(((datum-calc)/sigma)**2)
840        err = chisq/len(datum)
841        if report:
842            print ' Iteration: %d, chisq: %.3g, sum(shift^2): %.3g'%(nIter,chisq,np.sum(shift**2))
843    return chisq,Xw,calc
844
845###############################################################################
846#### SASD Utilities
847###############################################################################
848
849def SetScale(Data,refData):
850    Profile,Limits,Sample = Data
851    refProfile,refLimits,refSample = refData
852    x,y = Profile[:2]
853    rx,ry = refProfile[:2]
854    Beg = np.max([rx[0],x[0],Limits[1][0],refLimits[1][0]])
855    Fin = np.min([rx[-1],x[-1],Limits[1][1],refLimits[1][1]])
856    iBeg = np.searchsorted(x,Beg)
857    iFin = np.searchsorted(x,Fin)
858    sum = np.sum(y[iBeg:iFin])
859    refsum = np.sum(np.interp(x[iBeg:iFin],rx,ry,0,0))
860    Sample['Scale'][0] = refSample['Scale'][0]*refsum/sum
861   
862def Bestimate(G,Rg,P):
863    return (G*P/Rg**P)*np.exp(scsp.gammaln(P/2))
864   
865###############################################################################
866#### Size distribution
867###############################################################################
868
869def SizeDistribution(Profile,ProfDict,Limits,Substances,Sample,data):
870    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
871        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
872        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
873        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
874        'Cylinder diam':[CylinderDFF,CylinderDVol]}
875    Shape = data['Size']['Shape'][0]
876    Parms = data['Size']['Shape'][1:]
877    if data['Size']['logBins']:
878        Bins = np.logspace(np.log10(data['Size']['MinDiam']),np.log10(data['Size']['MaxDiam']),
879            data['Size']['Nbins']+1,True)/2.        #make radii
880    else:
881        Bins = np.linspace(data['Size']['MinDiam'],data['Size']['MaxDiam'],
882            data['Size']['Nbins']+1,True)/2.        #make radii
883    Dbins = np.diff(Bins)
884    Bins = Bins[:-1]+Dbins/2.
885    Contrast = Sample['Contrast'][1]
886    Scale = Sample['Scale'][0]
887    Sky = 10**data['Size']['MaxEnt']['Sky']
888    BinsBack = np.ones_like(Bins)*Sky*Scale/Contrast
889    Back = data['Back']
890    Q,Io,wt,Ic,Ib = Profile[:5]
891    Qmin = Limits[1][0]
892    Qmax = Limits[1][1]
893    wtFactor = ProfDict['wtFactor']
894    Ibeg = np.searchsorted(Q,Qmin)
895    Ifin = np.searchsorted(Q,Qmax)
896    BinMag = np.zeros_like(Bins)
897    Ic[:] = 0.
898    Gmat = G_matrix(Q[Ibeg:Ifin],Bins,Contrast,shapes[Shape][0],shapes[Shape][1],args=Parms)
899    if 'MaxEnt' == data['Size']['Method']:
900        chisq,BinMag,Ic[Ibeg:Ifin] = MaxEnt_SB(Scale*Io[Ibeg:Ifin]-Back[0],
901            Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),Gmat,BinsBack,
902            data['Size']['MaxEnt']['Niter'],report=True)
903    elif 'IPG' == data['Size']['Method']:
904        chisq,BinMag,Ic[Ibeg:Ifin] = IPG(Scale*Io[Ibeg:Ifin]-Back[0],Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),
905            Gmat,Bins,Dbins,data['Size']['IPG']['Niter'],Q[Ibeg:Ifin],approach=0.8,
906            Power=data['Size']['IPG']['Power'],report=True)
907    Ib[:] = Back[0]
908    Ic[Ibeg:Ifin] += Back[0]
909    print ' Final chi^2: %.3f'%(chisq)
910    Vols = shapes[Shape][1](Bins,Parms)
911    data['Size']['Distribution'] = [Bins,Dbins,BinMag/(2.*Dbins)]
912       
913################################################################################
914#### Modelling
915################################################################################
916
917def ModelFit(Profile,ProfDict,Limits,Substances,Sample,Model):
918    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
919        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
920        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
921        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
922        'Unified tube':[UniTubeFF,UniTubeVol],'Cylinder diam':[CylinderDFF,CylinderDVol]}
923
924    def GetModelParms():
925        parmOrder = ['Volume','Radius','Mean','StdDev','MinSize','G','Rg','B','P','Cutoff',
926            'PkInt','PkPos','PkSig','PkGam',]
927        parmDict = {'Scale':Sample['Scale'][0]}
928        varyList = []
929        values = []
930        levelTypes = []
931        Back = Model['Back']
932        if Back[1]:
933            varyList += ['Back',]
934            values.append(Back[0])
935        parmDict['Back'] = Back[0]
936        partData = Model['Particle']
937        parmDict['Matrix density'] = Substances['Substances'][partData['Matrix']['Name']].get('XAnom density',0.0)
938        for i,level in enumerate(partData['Levels']):
939            cid = str(i)+':'
940            controls = level['Controls']
941            Type = controls['DistType']
942            levelTypes.append(Type)
943            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
944                if Type not in ['Monodosperse',]:
945                    parmDict[cid+'NumPoints'] = controls['NumPoints']
946                    parmDict[cid+'Cutoff'] = controls['Cutoff']
947                parmDict[cid+'FormFact'] = shapes[controls['FormFact']][0]
948                parmDict[cid+'FFVolume'] = shapes[controls['FormFact']][1]
949                parmDict[cid+'XAnom density'] = Substances['Substances'][controls['Material']].get('XAnom density',0.0)
950                for item in ['Aspect ratio','Length','Thickness','Diameter',]:
951                    if item in controls['FFargs']:
952                        parmDict[cid+item] = controls['FFargs'][item][0]
953                        if controls['FFargs'][item][1]:
954                            varyList.append(cid+item)
955                            values.append(controls['FFargs'][item][0])
956            distDict = controls['DistType']
957            for item in parmOrder:
958                if item in level[distDict]:
959                    parmDict[cid+item] = level[distDict][item][0]
960                    if level[distDict][item][1]:
961                        values.append(level[distDict][item][0])
962                        varyList.append(cid+item)
963        return levelTypes,parmDict,varyList,values
964       
965    def SetModelParms():
966        print ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
967        if 'Back' in varyList:
968            Model['Back'][0] = parmDict['Back']
969            print '  %15s %15.4f esd: %15.4g'%('Background:',parmDict['Back'],sigDict['Back'])
970        partData = Model['Particle']
971        for i,level in enumerate(partData['Levels']):
972            Type = level['Controls']['DistType']
973            print ' Component %d: Type: %s:'%(i,Type)
974            cid = str(i)+':'
975            controls = level['Controls']
976            Type = controls['DistType']
977            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
978                for item in ['Aspect ratio','Length','Thickness','Diameter',]:
979                    if cid+item in varyList:
980                        controls['FFargs'][item][0] = parmDict[cid+item]
981                        print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
982            distDict = controls['DistType']
983            for item in level[distDict]:
984                if cid+item in varyList:
985                    level[distDict][item][0] = parmDict[cid+item]
986                    print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
987                   
988    def calcSASD(values,Q,Io,wt,Ifb,levelTypes,parmDict,varyList):
989        parmDict.update(zip(varyList,values))
990        M = np.sqrt(wt)*(getSASD(Q,levelTypes,parmDict)+Ifb-parmDict['Scale']*Io)
991        return M
992       
993    def getSASD(Q,levelTypes,parmDict):
994        Ic = np.zeros_like(Q)
995        rhoMat = parmDict['Matrix density']
996        for i,Type in enumerate(levelTypes):
997            cid = str(i)+':'
998            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm']:
999                FFfxn = parmDict[cid+'FormFact']
1000                Volfxn = parmDict[cid+'FFVolume']
1001                FFargs = []
1002                for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter',]:
1003                    if item in parmDict: 
1004                        FFargs.append(parmDict[item])
1005                distDict = {}
1006                for item in [cid+'Volume',cid+'Mean',cid+'StdDev',cid+'MinSize',]:
1007                    if item in parmDict:
1008                        distDict[item.split(':')[1]] = parmDict[item]
1009                rho = parmDict[cid+'XAnom density']
1010                contrast = rho**2-rhoMat**2
1011                rBins,dBins,dist = MakeDiamDist(Type,parmDict[cid+'NumPoints'],parmDict[cid+'Cutoff'],distDict)
1012                Gmat = G_matrix(Q,rBins,contrast,FFfxn,Volfxn,FFargs).T
1013                dist *= parmDict[cid+'Volume']
1014                Ic += np.dot(Gmat,dist)
1015            elif 'Unified' in Type:
1016                Rg,G,B,P,Rgco = parmDict[cid+'Rg'],parmDict[cid+'G'],parmDict[cid+'B'], \
1017                    parmDict[cid+'P'],parmDict[cid+'Cutoff']
1018                Qstar = Q/(scsp.erf(Q*Rg/np.sqrt(6)))**3
1019                Guin = G*np.exp(-(Q*Rg)**2/3)
1020                Porod = (B/Qstar**P)*np.exp(-(Q*Rgco)**2/3)
1021                Ic += Guin+Porod
1022            elif 'Porod' in Type:
1023                B,P,Rgco = parmDict[cid+'B'],parmDict[cid+'P'],parmDict[cid+'Cutoff']
1024                Porod = (B/Q**P)*np.exp(-(Q*Rgco)**2/3)
1025                Ic += Porod
1026            elif 'Mono' in Type:
1027                FFfxn = parmDict[cid+'FormFact']
1028                Volfxn = parmDict[cid+'FFVolume']
1029                FFargs = []
1030                for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter',]:
1031                    if item in parmDict: 
1032                        FFargs.append(parmDict[item])
1033                rho = parmDict[cid+'XAnom density']
1034                contrast = rho**2-rhoMat**2
1035                R = parmDict[cid+'Radius']
1036                Gmat = G_matrix(Q,R,contrast,FFfxn,Volfxn,FFargs)             
1037                Ic += Gmat[0]*parmDict[cid+'Volume']
1038            elif 'Bragg' in distFxn:
1039                Ic += parmDict[cid+'PkInt']*G2pwd.getPsVoigt(parmDict[cid+'PkPos'],
1040                    parmDict[cid+'PkSig'],parmDict[cid+'PkGam'],Q)
1041        Ic += parmDict['Back']  #/parmDict['Scale']
1042        return Ic
1043       
1044    Q,Io,wt,Ic,Ib,Ifb = Profile[:6]
1045    Qmin = Limits[1][0]
1046    Qmax = Limits[1][1]
1047    wtFactor = ProfDict['wtFactor']
1048    Ibeg = np.searchsorted(Q,Qmin)
1049    Ifin = np.searchsorted(Q,Qmax)
1050    Ic[:] = 0
1051    levelTypes,parmDict,varyList,values = GetModelParms()
1052    if varyList:
1053        result = so.leastsq(calcSASD,values,full_output=True,epsfcn=1.e-8,   #ftol=Ftol,
1054            args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Ifb[Ibeg:Ifin],levelTypes,parmDict,varyList))
1055        ncyc = int(result[2]['nfev']/2)
1056        parmDict.update(zip(varyList,result[0]))
1057        chisq = np.sum(result[2]['fvec']**2)
1058        ncalc = result[2]['nfev']
1059        covM = result[1]
1060    else:   #nothing varied
1061        M = calcSASD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Ifb[Ibeg:Ifin],levelTypes,parmDict,varyList)
1062        chisq = np.sum(M**2)
1063        ncalc = 0
1064        covM = []
1065        sig = []
1066        sigDict = {}
1067        result = []
1068    Rvals = {}
1069    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
1070    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
1071    Ic[Ibeg:Ifin] = getSASD(Q[Ibeg:Ifin],levelTypes,parmDict)
1072    try:
1073        if len(covM):
1074            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
1075            sigDict = dict(zip(varyList,sig))
1076        print ' Results of small angle data modelling fit:'
1077        print 'Number of function calls:',ncalc,' Number of observations: ',Ifin-Ibeg,' Number of parameters: ',len(varyList)
1078        print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
1079        SetModelParms()
1080        covMatrix = covM*Rvals['GOF']
1081        return True,result,varyList,sig,Rvals,covMatrix
1082    except ValueError:
1083        return False,0,0,0,0,0
1084   
1085def ModelFxn(Profile,ProfDict,Limits,Substances,Sample,sasdData):
1086    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
1087        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
1088        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
1089        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
1090        'Unified tube':[UniTubeFF,UniTubeVol],'Cylinder diam':[CylinderDFF,CylinderDVol]}
1091#    pdb.set_trace()
1092    partData = sasdData['Particle']
1093    rhoMat = Substances['Substances'][partData['Matrix']['Name']].get('XAnom density',0.0)
1094    matFrac = partData['Matrix']['VolFrac'] #[value,flag]       
1095    Scale = Sample['Scale']     #[value,flag]
1096    Back = sasdData['Back']
1097    Q,Io,wt,Ic,Ib,Ifb = Profile[:6]
1098    Qmin = Limits[1][0]
1099    Qmax = Limits[1][1]
1100    wtFactor = ProfDict['wtFactor']
1101    Ibeg = np.searchsorted(Q,Qmin)
1102    Ifin = np.searchsorted(Q,Qmax)
1103    Ib[:] = Back[0]
1104    Ic[:] = 0
1105    Rbins = []
1106    Dist = []
1107    for level in partData['Levels']:
1108        controls = level['Controls']
1109        distFxn = controls['DistType']
1110        if distFxn in ['LogNormal','Gaussian','LSW','Schulz-Zimm']:
1111            FFfxn = shapes[controls['FormFact']][0]
1112            Volfxn = shapes[controls['FormFact']][1]
1113            FFargs = []
1114            for item in ['Aspect ratio','Length','Thickness','Diameter',]:
1115                if item in controls['FFargs']: 
1116                    FFargs.append(controls['FFargs'][item][0])
1117            rho = Substances['Substances'][level['Controls']['Material']].get('XAnom density',0.0)
1118            contrast = rho**2-rhoMat**2
1119            parmDict = level[controls['DistType']]
1120            distDict = {}
1121            for item in parmDict:
1122                distDict[item] = parmDict[item][0]
1123            rBins,dBins,dist = MakeDiamDist(controls['DistType'],controls['NumPoints'],controls['Cutoff'],distDict)
1124            Gmat = G_matrix(Q[Ibeg:Ifin],rBins,contrast,FFfxn,Volfxn,FFargs).T
1125            dist *= level[distFxn]['Volume'][0]
1126            Ic[Ibeg:Ifin] += np.dot(Gmat,dist)
1127            Rbins.append(rBins)
1128            Dist.append(dist/(4.*dBins))
1129        elif 'Unified' in distFxn:
1130            parmDict = level[controls['DistType']]
1131            Rg,G,B,P,Rgco = parmDict['Rg'][0],parmDict['G'][0],parmDict['B'][0],    \
1132                parmDict['P'][0],parmDict['Cutoff'][0]
1133            Qstar = Q[Ibeg:Ifin]/(scsp.erf(Q[Ibeg:Ifin]*Rg/np.sqrt(6)))**3
1134            Guin = G*np.exp(-(Q[Ibeg:Ifin]*Rg)**2/3)
1135            Porod = (B/Qstar**P)*np.exp(-(Q[Ibeg:Ifin]*Rgco)**2/3)
1136            Ic[Ibeg:Ifin] += Guin+Porod
1137            Rbins.append([])
1138            Dist.append([])
1139        elif 'Porod' in distFxn:
1140            parmDict = level[controls['DistType']]
1141            B,P,Rgco = parmDict['B'][0],parmDict['P'][0],parmDict['Cutoff'][0]
1142            Porod = (B/Q[Ibeg:Ifin]**P)*np.exp(-(Q[Ibeg:Ifin]*Rgco)**2/3)
1143            Ic[Ibeg:Ifin] += Porod
1144            Rbins.append([])
1145            Dist.append([])
1146        elif 'Mono' in distFxn:
1147            FFfxn = shapes[controls['FormFact']][0]
1148            Volfxn = shapes[controls['FormFact']][1]
1149            FFargs = []
1150            for item in ['Aspect ratio','Length','Thickness','Diameter',]:
1151                if item in controls['FFargs']: 
1152                    FFargs.append(controls['FFargs'][item][0])
1153            rho = Substances['Substances'][level['Controls']['Material']].get('XAnom density',0.0)
1154            contrast = rho**2-rhoMat**2
1155            R = level[controls['DistType']]['Radius'][0]
1156            Gmat = G_matrix(Q[Ibeg:Ifin],R,contrast,FFfxn,Volfxn,FFargs)             
1157            Ic[Ibeg:Ifin] += Gmat[0]*level[distFxn]['Volume'][0]
1158            Rbins.append([])
1159            Dist.append([])
1160        elif 'Bragg' in distFxn:
1161            parmDict = level[controls['DistType']]
1162            Ic[Ibeg:Ifin] += parmDict['PkInt'][0]*G2pwd.getPsVoigt(parmDict['PkPos'][0],
1163                parmDict['PkSig'][0],parmDict['PkGam'][0],Q[Ibeg:Ifin])
1164            Rbins.append([])
1165            Dist.append([])
1166    Ic[Ibeg:Ifin] += Back[0]
1167    sasdData['Size Calc'] = [Rbins,Dist]
1168   
1169def MakeDiamDist(DistName,nPoints,cutoff,distDict):
1170   
1171    if 'LogNormal' in DistName:
1172        distFxn = 'LogNormalDist'
1173        cumeFxn = 'LogNormalCume'
1174        pos = distDict['MinSize']
1175        args = [distDict['Mean'],distDict['StdDev']]
1176        step = 4.*np.sqrt(np.exp(distDict['StdDev']**2)*(np.exp(distDict['StdDev']**2)-1.))
1177        mode = distDict['MinSize']+distDict['Mean']/np.exp(distDict['StdDev']**2)
1178        minX = 1. #pos
1179        maxX = 1.e4 #np.min([mode+nPoints*step*40,1.e5])
1180    elif 'Gauss' in DistName:
1181        distFxn = 'GaussDist'
1182        cumeFxn = 'GaussCume'
1183        pos = distDict['Mean']
1184        args = [distDict['StdDev']]
1185        step = 0.02*distDict['StdDev']
1186        mode = distDict['Mean']
1187        minX = np.max([mode-4.*distDict['StdDev'],1.])
1188        maxX = np.min([mode+4.*distDict['StdDev'],1.e5])
1189    elif 'LSW' in DistName:
1190        distFxn = 'LSWDist'
1191        cumeFxn = 'LSWCume'
1192        pos = distDict['Mean']
1193        args = []
1194        minX,maxX = [0.,2.*pos]
1195    elif 'Schulz' in DistName:
1196        distFxn = 'SchulzZimmDist'
1197        cumeFxn = 'SchulzZimmCume'
1198        pos = distDict['Mean']
1199        args = [distDict['StdDev']]
1200        minX = np.max([1.,pos-4.*distDict['StdDev']])
1201        maxX = np.min([pos+4.*distDict['StdDev'],1.e5])
1202    nP = 500
1203    Diam = np.logspace(0.,5.,nP,True)
1204    TCW = eval(cumeFxn+'(Diam,pos,args)')
1205    CumeTgts = np.linspace(cutoff,(1.-cutoff),nPoints+1,True)
1206    rBins = np.interp(CumeTgts,TCW,Diam,0,0)
1207    dBins = np.diff(rBins)
1208    rBins = rBins[:-1]+dBins/2.
1209    return rBins,dBins,eval(distFxn+'(rBins,pos,args)')
1210
1211################################################################################
1212#### MaxEnt testing stuff
1213################################################################################
1214
1215def print_vec(text, a):
1216    '''print the contents of a vector to the console'''
1217    n = a.shape[0]
1218    print "%s[ = (" % text,
1219    for i in range(n):
1220        s = " %g, " % a[i]
1221        print s,
1222    print ")"
1223
1224def print_arr(text, a):
1225    '''print the contents of an array to the console'''
1226    n, m = a.shape
1227    print "%s[][] = (" % text
1228    for i in range(n):
1229        print " (",
1230        for j in range(m):
1231            print " %g, " % a[i][j],
1232        print "),"
1233    print ")"
1234
1235def test_MaxEnt_SB(report=True):
1236    def readTextData(filename):
1237        '''return q, I, dI from a 3-column text file'''
1238        if not os.path.exists(filename):
1239            raise Exception("file not found: " + filename)
1240        buf = [line.split() for line in open(filename, 'r').readlines()]
1241        M = len(buf)
1242        buf = zip(*buf)         # transpose rows and columns
1243        q  = np.array(buf[0], dtype=np.float64)
1244        I  = np.array(buf[1], dtype=np.float64)
1245        dI = np.array(buf[2], dtype=np.float64)
1246        return q, I, dI
1247    print "MaxEnt_SB: "
1248    test_data_file = os.path.join( 'testinp', 'test.sas')
1249    rhosq = 100     # scattering contrast, 10^20 1/cm^-4
1250    bkg   = 0.1     #   I = I - bkg
1251    dMin, dMax, nRadii = 25, 9000, 40
1252    defaultDistLevel = 1.0e-6
1253    IterMax = 40
1254    errFac = 1.05
1255   
1256    r    = np.logspace(math.log10(dMin), math.log10(dMax), nRadii)/2
1257    dr   = r * (r[1]/r[0] - 1)          # step size
1258    f_dr = np.ndarray((nRadii)) * 0  # volume fraction histogram
1259    b    = np.ndarray((nRadii)) * 0 + defaultDistLevel  # MaxEnt "sky background"
1260   
1261    qVec, I, dI = readTextData(test_data_file)
1262    G = G_matrix(qVec,r,rhosq,SphereFF,SphereVol,args=())
1263   
1264    chisq,f_dr,Ic = MaxEnt_SB(I - bkg, dI*errFac, b, IterMax, G, report=report)
1265    if f_dr is None:
1266        print "no solution"
1267        return
1268   
1269    print "solution reached"
1270    for a,b,c in zip(r.tolist(), dr.tolist(), f_dr.tolist()):
1271        print '%10.4f %10.4f %12.4g'%(a,b,c)
1272
1273def tests():
1274    test_MaxEnt_SB(report=True)
1275
1276if __name__ == '__main__':
1277    tests()
1278
Note: See TracBrowser for help on using the repository browser.