Changeset 1302


Ignore:
Timestamp:
Apr 25, 2014 12:56:35 PM (10 years ago)
Author:
vondreele
Message:

merge messup fixed!

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIsasd.py

    r1301 r1302  
    1 <<<<<<< .mine
    21#/usr/bin/env python
    32# -*- coding: utf-8 -*-
    … …  
    12781277    tests()
    12791278
    1280 =======
    1281 #/usr/bin/env python
    1282 # -*- coding: utf-8 -*-
    1283 '''
    1284 *GSASII small angle calculation module*
    1285 ========================================
    1286 
    1287 '''
    1288 ########### SVN repository information ###################
    1289 # $Date$
    1290 # $Author$
    1291 # $Revision$
    1292 # $URL$
    1293 # $Id$
    1294 ########### SVN repository information ###################
    1295 import os
    1296 import sys
    1297 import math
    1298 import time
    1299 
    1300 import numpy as np
    1301 import scipy as sp
    1302 import numpy.linalg as nl
    1303 from numpy.fft import ifft, fft, fftshift
    1304 import scipy.special as scsp
    1305 import scipy.interpolate as si
    1306 import scipy.stats as st
    1307 import scipy.optimize as so
    1308 #import pdb
    1309 
    1310 import GSASIIpath
    1311 GSASIIpath.SetVersionNumber("$Revision$")
    1312 import GSASIIlattice as G2lat
    1313 import GSASIIspc as G2spc
    1314 import GSASIIElem as G2elem
    1315 import GSASIIgrid as G2gd
    1316 import GSASIIIO as G2IO
    1317 import GSASIImath as G2mth
    1318 import GSASIIpwd as G2pwd
    1319 
    1320 # trig functions in degrees
    1321 sind = lambda x: math.sin(x*math.pi/180.)
    1322 asind = lambda x: 180.*math.asin(x)/math.pi
    1323 tand = lambda x: math.tan(x*math.pi/180.)
    1324 atand = lambda x: 180.*math.atan(x)/math.pi
    1325 atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
    1326 cosd = lambda x: math.cos(x*math.pi/180.)
    1327 acosd = lambda x: 180.*math.acos(x)/math.pi
    1328 rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
    1329 #numpy versions
    1330 npsind = lambda x: np.sin(x*np.pi/180.)
    1331 npasind = lambda x: 180.*np.arcsin(x)/math.pi
    1332 npcosd = lambda x: np.cos(x*math.pi/180.)
    1333 npacosd = lambda x: 180.*np.arccos(x)/math.pi
    1334 nptand = lambda x: np.tan(x*math.pi/180.)
    1335 npatand = lambda x: 180.*np.arctan(x)/np.pi
    1336 npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
    1337 npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave
    1338 npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
    1339    
    1340 ###############################################################################
    1341 #### Particle form factors
    1342 ###############################################################################
    1343 
    1344 def SphereFF(Q,R,args=()):
    1345     ''' Compute hard sphere form factor - can use numpy arrays
    1346     param float Q: Q value array (usually in A-1)
    1347     param float R: sphere radius (Usually in A - must match Q-1 units)
    1348     param array args: ignored
    1349     returns float: form factors as array as needed
    1350     '''
    1351     QR = Q[:,np.newaxis]*R
    1352     return (3./(QR**3))*(np.sin(QR)-(QR*np.cos(QR)))
    1353    
    1354 def SpheroidFF(Q,R,args):
    1355     ''' Compute form factor of cylindrically symmetric ellipsoid (spheroid)
    1356     - can use numpy arrays for R & AR; will return corresponding numpy array
    1357     param float Q : Q value array (usually in A-1)
    1358     param float R: radius along 2 axes of spheroid
    1359     param array args: [float AR]: aspect ratio so 3rd axis = R*AR
    1360     returns float: form factors as array as needed
    1361     '''
    1362     NP = 50
    1363     AR = args[0]
    1364     if 0.99 < AR < 1.01:
    1365         return SphereFF(Q,R,0)
    1366     else:
    1367         cth = np.linspace(0,1.,NP)
    1368         Rct = R[:,np.newaxis]*np.sqrt(1.+(AR**2-1.)*cth**2)
    1369         return np.sqrt(np.sum(SphereFF(Q[:,np.newaxis],Rct,0)**2,axis=2)/NP)
    1370            
    1371 def CylinderFF(Q,R,args):
    1372     ''' Compute form factor for cylinders - can use numpy arrays
    1373     param float Q: Q value array (A-1)
    1374     param float R: cylinder radius (A)
    1375     param array args: [float L]: cylinder length (A)
    1376     returns float: form factor
    1377     '''
    1378     L = args[0]
    1379     NP = 200
    1380     alp = np.linspace(0,np.pi/2.,NP)
    1381     QL = Q[:,np.newaxis]*L
    1382     LBessArg = 0.5*QL[:,:,np.newaxis]**np.cos(alp)
    1383     LBess = np.where(LBessArg<1.e-6,1.,np.sin(LBessArg)/LBessArg)
    1384     QR = Q[:,np.newaxis]*R
    1385     SBessArg = QR[:,:,np.newaxis]*np.sin(alp)
    1386     SBess = np.where(SBessArg<1.e-6,0.5,scsp.jv(1,SBessArg)/SBessArg)
    1387     LSBess = LBess*SBess
    1388     return np.sqrt(2.*np.pi*np.sum(np.sin(alp)*LSBess**2,axis=2)/NP)
    1389    
    1390 def CylinderDFF(Q,L,args):
    1391     ''' Compute form factor for cylinders - can use numpy arrays
    1392     param float Q: Q value array (A-1)
    1393     param float L: cylinder half length (A)
    1394     param array args: [float R]: cylinder radius (A)
    1395     returns float: form factor
    1396     '''
    1397     R = args[0]
    1398     return CylinderFF(Q,R,args=[2.*L])   
    1399    
    1400 def CylinderARFF(Q,R,args):
    1401     ''' Compute form factor for cylinders - can use numpy arrays
    1402     param float Q: Q value array (A-1)
    1403     param float R: cylinder radius (A)
    1404     param array args: [float AR]: cylinder aspect ratio = L/D = L/2R
    1405     returns float: form factor
    1406     '''
    1407     AR = args[0]
    1408     return CylinderFF(Q,R,args=[2.*R*AR])   
    1409    
    1410 def UniSphereFF(Q,R,args=0):
    1411     ''' Compute form factor for unified sphere - can use numpy arrays
    1412     param float Q: Q value array (A-1)
    1413     param float R: cylinder radius (A)
    1414     param array args: ignored
    1415     returns float: form factor
    1416     '''
    1417     Rg = np.sqrt(3./5.)*R
    1418     B = np.pi*1.62/(Rg**4)    #are we missing *np.pi? 1.62 = 6*(3/5)**2/(4/3) sense?
    1419     QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg/np.sqrt(6)))**3
    1420     return np.sqrt(np.exp((-Q[:,np.newaxis]**2*Rg**2)/3.)+(B/QstV**4))
    1421    
    1422 def UniRodFF(Q,R,args):
    1423     ''' Compute form factor for unified rod - can use numpy arrays
    1424     param float Q: Q value array (A-1)
    1425     param float R: cylinder radius (A)
    1426     param array args: [float R]: cylinder radius (A)
    1427     returns float: form factor
    1428     '''
    1429     L = args[0]
    1430     Rg2 = np.sqrt(R**2/2+L**2/12)
    1431     B2 = np.pi/L
    1432     Rg1 = np.sqrt(3.)*R/2.
    1433     G1 = (2./3.)*R/L
    1434     B1 = 4.*(L+R)/(R**3*L**2)
    1435     QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg2/np.sqrt(6)))**3
    1436     FF = np.exp(-Q[:,np.newaxis]**2*Rg2**2/3.)+(B2/QstV)*np.exp(-Rg1**2*Q[:,np.newaxis]**2/3.)
    1437     QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg1/np.sqrt(6)))**3
    1438     FF += G1*np.exp(-Q[:,np.newaxis]**2*Rg1**2/3.)+(B1/QstV**4)
    1439     return np.sqrt(FF)
    1440    
    1441 def UniRodARFF(Q,R,args):
    1442     ''' Compute form factor for unified rod of fixed aspect ratio - can use numpy arrays
    1443     param float Q: Q value array (A-1)
    1444     param float R: cylinder radius (A)
    1445     param array args: [float AR]: cylinder aspect ratio = L/D = L/2R
    1446     returns float: form factor
    1447     '''
    1448     AR = args[0]
    1449     return UniRodFF(Q,R,args=[2.*AR*R,])
    1450    
    1451 def UniDiskFF(Q,R,args):
    1452     ''' Compute form factor for unified disk - can use numpy arrays
    1453     param float Q: Q value array (A-1)
    1454     param float R: cylinder radius (A)
    1455     param array args: [float T]: disk thickness (A)
    1456     returns float: form factor
    1457     '''
    1458     T = args[0]
    1459     Rg2 = np.sqrt(R**2/2.+T**2/12.)
    1460     B2 = 2./R**2
    1461     Rg1 = np.sqrt(3.)*T/2.
    1462     RgC2 = 1.1*Rg1
    1463     G1 = (2./3.)*(T/R)**2
    1464     B1 = 4.*(T+R)/(R**3*T**2)
    1465     QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg2/np.sqrt(6)))**3
    1466     FF = np.exp(-Q[:,np.newaxis]**2*Rg2**2/3.)+(B2/QstV**2)*np.exp(-RgC2**2*Q[:,np.newaxis]**2/3.)
    1467     QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg1/np.sqrt(6)))**3
    1468     FF += G1*np.exp(-Q[:,np.newaxis]**2*Rg1**2/3.)+(B1/QstV**4)
    1469     return np.sqrt(FF)
    1470    
    1471 def UniTubeFF(Q,R,args):
    1472     ''' Compute form factor for unified tube - can use numpy arrays
    1473     assumes that core of tube is same as the matrix/solvent so contrast
    1474     is from tube wall vs matrix
    1475     param float Q: Q value array (A-1)
    1476     param float R: cylinder radius (A)
    1477     param array args: [float L,T]: tube length & wall thickness(A)
    1478     returns float: form factor
    1479     '''
    1480     L,T = args[:2]
    1481     Ri = R-T
    1482     DR2 = R**2-Ri**2
    1483     Vt = np.pi*DR2*L
    1484     Rg3 = np.sqrt(DR2/2.+L**2/12.)
    1485     B1 = 4.*np.pi**2*(DR2+L*(R+Ri))/Vt**2
    1486     B2 = np.pi**2*T/Vt
    1487     B3 = np.pi/L
    1488     QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg3/np.sqrt(6)))**3
    1489     FF = np.exp(-Q[:,np.newaxis]**2*Rg3**2/3.)+(B3/QstV)*np.exp(-Q[:,np.newaxis]**2*R**2/3.)
    1490     QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*R/np.sqrt(6)))**3
    1491     FF += (B2/QstV**2)*np.exp(-Q[:,np.newaxis]**2*T**2/3.)
    1492     QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*T/np.sqrt(6)))**3
    1493     FF += B1/QstV**4
    1494     return np.sqrt(FF)
    1495 
    1496 ###############################################################################
    1497 #### Particle volumes
    1498 ###############################################################################
    1499 
    1500 def SphereVol(R,args=()):
    1501     ''' Compute volume of sphere
    1502     - numpy array friendly
    1503     param float R: sphere radius
    1504     param array args: ignored
    1505     returns float: volume
    1506     '''
    1507     return (4./3.)*np.pi*R**3
    1508 
    1509 def SpheroidVol(R,args):
    1510     ''' Compute volume of cylindrically symmetric ellipsoid (spheroid)
    1511     - numpy array friendly
    1512     param float R: radius along 2 axes of spheroid
    1513     param array args: [float AR]: aspect ratio so radius of 3rd axis = R*AR
    1514     returns float: volume
    1515     '''
    1516     AR = args[0]
    1517     return AR*SphereVol(R)
    1518    
    1519 def CylinderVol(R,args):
    1520     ''' Compute cylinder volume for radius & length
    1521     - numpy array friendly
    1522     param float R: diameter (A)
    1523     param array args: [float L]: length (A)
    1524     returns float:volume (A^3)
    1525     '''
    1526     L = args[0]
    1527     return np.pi*L*R**2
    1528    
    1529 def CylinderDVol(L,args):
    1530     ''' Compute cylinder volume for length & diameter
    1531     - numpy array friendly
    1532     param float: L half length (A)
    1533     param array args: [float D]: diameter (A)
    1534     returns float:volume (A^3)
    1535     '''
    1536     D = args[0]
    1537     return CylinderVol(D/2.,[2.*L,])
    1538    
    1539 def CylinderARVol(R,args):
    1540     ''' Compute cylinder volume for radius & aspect ratio = L/D
    1541     - numpy array friendly
    1542     param float: R radius (A)
    1543     param array args: [float AR]: =L/D=L/2R aspect ratio
    1544     returns float:volume
    1545     '''
    1546     AR = args[0]
    1547     return CylinderVol(R,[2.*R*AR,])
    1548    
    1549 def UniSphereVol(R,args=()):
    1550     ''' Compute volume of sphere
    1551     - numpy array friendly
    1552     param float R: sphere radius
    1553     param array args: ignored
    1554     returns float: volume
    1555     '''
    1556     return SphereVol(R)
    1557    
    1558 def UniRodVol(R,args):
    1559     ''' Compute cylinder volume for radius & length
    1560     - numpy array friendly
    1561     param float R: diameter (A)
    1562     param array args: [float L]: length (A)
    1563     returns float:volume (A^3)
    1564     '''
    1565     L = args[0]
    1566     return CylinderVol(R,[L,])
    1567    
    1568 def UniRodARVol(R,args):
    1569     ''' Compute rod volume for radius & aspect ratio
    1570     - numpy array friendly
    1571     param float R: diameter (A)
    1572     param array args: [float AR]: =L/D=L/2R aspect ratio
    1573     returns float:volume (A^3)
    1574     '''
    1575     AR = args[0]
    1576     return CylinderARVol(R,[AR,])
    1577    
    1578 def UniDiskVol(R,args):
    1579     ''' Compute disk volume for radius & thickness
    1580     - numpy array friendly
    1581     param float R: diameter (A)
    1582     param array args: [float T]: thickness
    1583     returns float:volume (A^3)
    1584     '''
    1585     T = args[0]
    1586     return CylinderVol(R,[T,])
    1587    
    1588 def UniTubeVol(R,args):
    1589     ''' Compute tube volume for radius, length & wall thickness
    1590     - numpy array friendly
    1591     param float R: diameter (A)
    1592     param array args: [float L,T]: tube length & wall thickness(A)
    1593     returns float: volume (A^3) of tube wall
    1594     '''
    1595     L,T = args[:2]
    1596     return CylinderVol(R,[L,])-CylinderVol(R-T,[L,])
    1597    
    1598 ################################################################################
    1599 #### Distribution functions & their cumulative fxns
    1600 ################################################################################
    1601 
    1602 def LogNormalDist(x,pos,args):
    1603     ''' Standard LogNormal distribution - numpy friendly on x axis
    1604     ref: http://www.itl.nist.gov/div898/handbook/index.htm 1.3.6.6.9
    1605     param float x: independent axis (can be numpy array)
    1606     param float pos: location of distribution
    1607     param float scale: width of distribution (m)
    1608     param float shape: shape - (sigma of log(LogNormal))
    1609     returns float: LogNormal distribution
    1610     '''
    1611     scale,shape = args
    1612     return np.exp(-np.log((x-pos)/scale)**2/(2.*shape**2))/(np.sqrt(2.*np.pi)*(x-pos)*shape)
    1613    
    1614 def GaussDist(x,pos,args):
    1615     ''' Standard Normal distribution - numpy friendly on x axis
    1616     param float x: independent axis (can be numpy array)
    1617     param float pos: location of distribution
    1618     param float scale: width of distribution (sigma)
    1619     param float shape: not used
    1620     returns float: Normal distribution
    1621     '''
    1622     scale = args[0]
    1623     return (1./(scale*np.sqrt(2.*np.pi)))*np.exp(-(x-pos)**2/(2.*scale**2))
    1624    
    1625 def LSWDist(x,pos,args=[]):
    1626     ''' Lifshitz-Slyozov-Wagner Ostwald ripening distribution - numpy friendly on x axis
    1627     ref:
    1628     param float x: independent axis (can be numpy array)
    1629     param float pos: location of distribution
    1630     param float scale: not used
    1631     param float shape: not used
    1632     returns float: LSW distribution   
    1633     '''
    1634     redX = x/pos
    1635     result = (81.*2**(-5/3.))*(redX**2*np.exp(-redX/(1.5-redX)))/((1.5-redX)**(11/3.)*(3.-redX)**(7/3.))
    1636    
    1637     return np.nan_to_num(result/pos)
    1638    
    1639 def SchulzZimmDist(x,pos,args):
    1640     ''' Schulz-Zimm macromolecule distribution - numpy friendly on x axis
    1641     ref: http://goldbook.iupac.org/S05502.html
    1642     param float x: independent axis (can be numpy array)
    1643     param float pos: location of distribution
    1644     param float scale: width of distribution (sigma)
    1645     param float shape: not used
    1646     returns float: Schulz-Zimm distribution
    1647     '''
    1648     scale = args[0]
    1649     b = (2.*pos/scale)**2
    1650     a = b/pos
    1651     if b<70:    #why bother?
    1652         return (a**(b+1.))*x**b*np.exp(-a*x)/scsp.gamma(b+1.)
    1653     else:
    1654         return np.exp((b+1.)*np.log(a)-scsp.gammaln(b+1.)+b*np.log(x)-(a*x))
    1655            
    1656 def LogNormalCume(x,pos,args):
    1657     ''' Standard LogNormal cumulative distribution - numpy friendly on x axis
    1658     ref: http://www.itl.nist.gov/div898/handbook/index.htm 1.3.6.6.9
    1659     param float x: independent axis (can be numpy array)
    1660     param float pos: location of distribution
    1661     param float scale: width of distribution (sigma)
    1662     param float shape: shape parameter
    1663     returns float: LogNormal cumulative distribution
    1664     '''
    1665     scale,shape = args
    1666     lredX = np.log((x-pos)/scale)
    1667     return (scsp.erf((lredX/shape)/np.sqrt(2.))+1.)/2.
    1668    
    1669 def GaussCume(x,pos,args):
    1670     ''' Standard Normal cumulative distribution - numpy friendly on x axis
    1671     param float x: independent axis (can be numpy array)
    1672     param float pos: location of distribution
    1673     param float scale: width of distribution (sigma)
    1674     param float shape: not used
    1675     returns float: Normal cumulative distribution
    1676     '''
    1677     scale = args[0]
    1678     redX = (x-pos)/scale
    1679     return (scsp.erf(redX/np.sqrt(2.))+1.)/2.
    1680    
    1681 def LSWCume(x,pos,args=[]):
    1682     ''' Lifshitz-Slyozov-Wagner Ostwald ripening cumulative distribution - numpy friendly on x axis
    1683     param float x: independent axis (can be numpy array)
    1684     param float pos: location of distribution
    1685     param float scale: not used
    1686     param float shape: not used
    1687     returns float: LSW cumulative distribution
    1688     '''
    1689     nP = 500
    1690     xMin,xMax = [0.,2.*pos]
    1691     X = np.linspace(xMin,xMax,nP,True)
    1692     fxn = LSWDist(X,pos)
    1693     mat = np.outer(np.ones(nP),fxn)
    1694     cume = np.sum(np.tril(mat),axis=1)/np.sum(fxn)
    1695     return np.interp(x,X,cume,0,1)
    1696    
    1697 def SchulzZimmCume(x,pos,args):
    1698     ''' Schulz-Zimm cumulative distribution - numpy friendly on x axis
    1699     param float x: independent axis (can be numpy array)
    1700     param float pos: location of distribution
    1701     param float scale: width of distribution (sigma)
    1702     param float shape: not used
    1703     returns float: Normal distribution
    1704     '''
    1705     scale = args[0]
    1706     nP = 500
    1707     xMin = np.max([0.,pos-4.*scale])
    1708     xMax = np.min([pos+4.*scale,1.e5])
    1709     X = np.linspace(xMin,xMax,nP,True)
    1710     fxn = LSWDist(X,pos)
    1711     mat = np.outer(np.ones(nP),fxn)
    1712     cume = np.sum(np.tril(mat),axis=1)/np.sum(fxn)
    1713     return np.interp(x,X,cume,0,1)
    1714    
    1715     return []
    1716    
    1717          
    1718 ################################################################################
    1719 ##### SB-MaxEnt
    1720 ################################################################################
    1721 
    1722 def G_matrix(q,r,contrast,FFfxn,Volfxn,args=()):
    1723     '''Calculates the response matrix :math:`G(Q,r)`
    1724    
    1725     :param float q: :math:`Q`
    1726     :param float r: :math:`r`
    1727     :param float contrast: :math:`|\\Delta\\rho|^2`, the scattering contrast
    1728     :param function FFfxn: form factor function FF(q,r,args)
    1729     :param function Volfxn: volume function Vol(r,args)
    1730     :returns float: G(Q,r)
    1731     '''
    1732     FF = FFfxn(q,r,args)
    1733     Vol = Volfxn(r,args)
    1734     return 1.e-4*(contrast*Vol*FF**2).T     #10^-20 vs 10^-24
    1735    
    1736 '''
    1737 sbmaxent
    1738 
    1739 Entropy maximization routine as described in the article
    1740 J Skilling and RK Bryan; MNRAS 211 (1984) 111 - 124.
    1741 ("MNRAS": "Monthly Notices of the Royal Astronomical Society")
    1742 
    1743 :license: Copyright (c) 2013, UChicago Argonne, LLC
    1744 :license: This file is distributed subject to a Software License Agreement found
    1745      in the file LICENSE that is included with this distribution.
    1746 
    1747 References:
    1748 
    1749 1. J Skilling and RK Bryan; MON NOT R ASTR SOC 211 (1984) 111 - 124.
    1750 2. JA Potton, GJ Daniell, and BD Rainford; Proc. Workshop
    1751    Neutron Scattering Data Analysis, Rutherford
    1752    Appleton Laboratory, UK, 1986; ed. MW Johnson,
    1753    IOP Conference Series 81 (1986) 81 - 86, Institute
    1754    of Physics, Bristol, UK.
    1755 3. ID Culverwell and GP Clarke; Ibid. 87 - 96.
    1756 4. JA Potton, GK Daniell, & BD Rainford,
    1757    J APPL CRYST 21 (1988) 663 - 668.
    1758 5. JA Potton, GJ Daniell, & BD Rainford,
    1759    J APPL CRYST 21 (1988) 891 - 897.
    1760 
    1761 '''
    1762 
    1763 class MaxEntException(Exception):
    1764     '''Any exception from this module'''
    1765     pass
    1766 
    1767 def MaxEnt_SB(datum, sigma, G, base, IterMax, image_to_data=None, data_to_image=None, report=False):
    1768     '''
    1769     do the complete Maximum Entropy algorithm of Skilling and Bryan
    1770    
    1771     :param float datum[]:
    1772     :param float sigma[]:
    1773     :param float[][] G: transformation matrix
    1774     :param float base[]:
    1775     :param int IterMax:
    1776     :param obj image_to_data: opus function (defaults to opus)
    1777     :param obj data_to_image: tropus function (defaults to tropus)
    1778    
    1779     :returns float[]: :math:`f(r) dr`
    1780     '''
    1781        
    1782     TEST_LIMIT        = 0.05                    # for convergence
    1783     CHI_SQR_LIMIT     = 0.01                    # maximum difference in ChiSqr for a solution
    1784     SEARCH_DIRECTIONS = 3                       # <10.  This code requires value = 3
    1785     RESET_STRAYS      = 1                       # was 0.001, correction of stray negative values
    1786     DISTANCE_LIMIT_FACTOR = 0.1                 # limitation on df to constrain runaways
    1787    
    1788     MAX_MOVE_LOOPS    = 500                     # for no solution in routine: move,
    1789     MOVE_PASSES       = 0.001                   # convergence test in routine: move
    1790 
    1791     def tropus (data, G):
    1792         '''
    1793         tropus: transform data-space -> solution-space:  [G] * data
    1794        
    1795         default definition, caller can use this definition or provide an alternative
    1796        
    1797         :param float[M] data: observations, ndarray of shape (M)
    1798         :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
    1799         :returns float[N]: calculated image, ndarray of shape (N)
    1800         '''
    1801         return G.dot(data)
    1802 
    1803     def opus (image, G):
    1804         '''
    1805         opus: transform solution-space -> data-space:  [G]^tr * image
    1806        
    1807         default definition, caller can use this definition or provide an alternative
    1808        
    1809         :param float[N] image: solution, ndarray of shape (N)
    1810         :param float[M][N] G: transformation matrix, ndarray of shape (M,N)
    1811         :returns float[M]: calculated data, ndarray of shape (M)
    1812         '''
    1813         return np.dot(G.T,image)    #G.transpose().dot(image)
    1814 
    1815     def Dist(s2, beta):
    1816         '''measure the distance of this possible solution'''
    1817         w = 0
    1818         n = beta.shape[0]
    1819         for k in range(n):
    1820             z = -sum(s2[k] * beta)
    1821             w += beta[k] * z
    1822         return w
    1823    
    1824     def ChiNow(ax, c1, c2, s1, s2):
    1825         '''
    1826         ChiNow
    1827        
    1828         :returns tuple: (ChiNow computation of ``w``, beta)
    1829         '''
    1830        
    1831         bx = 1 - ax
    1832         a =   bx * c2  -  ax * s2
    1833         b = -(bx * c1  -  ax * s1)
    1834    
    1835         beta = ChoSol(a, b)
    1836         w = 1.0
    1837         for k in range(SEARCH_DIRECTIONS):
    1838             w += beta[k] * (c1[k] + 0.5*sum(c2[k] * beta))
    1839         return w, beta
    1840    
    1841     def ChoSol(a, b):
    1842         '''
    1843         ChoSol: ? chop the solution vectors ?
    1844        
    1845         :returns: new vector beta
    1846         '''
    1847         n = b.shape[0]
    1848         fl = np.zeros((n,n))
    1849         bl = np.zeros_like(b)
    1850        
    1851         #print_arr("ChoSol: a", a)
    1852         #print_vec("ChoSol: b", b)
    1853    
    1854         if (a[0][0] <= 0):
    1855             msg = "ChoSol: a[0][0] = "
    1856             msg += str(a[0][0])
    1857             msg += '  Value must be positive'
    1858             raise MaxEntException(msg)
    1859    
    1860         # first, compute fl from a
    1861         # note fl is a lower triangular matrix
    1862         fl[0][0] = math.sqrt (a[0][0])
    1863         for i in (1, 2):
    1864             fl[i][0] = a[i][0] / fl[0][0]
    1865             for j in range(1, i+1):
    1866                 z = 0.0
    1867                 for k in range(j):
    1868                     z += fl[i][k] * fl[j][k]
    1869                     #print "ChoSol: %d %d %d  z = %lg" % ( i, j, k, z)
    1870                 z = a[i][j] - z
    1871                 if j == i:
    1872                     y = math.sqrt(z)
    1873                 else:
    1874                     y = z / fl[j][j]
    1875                 fl[i][j] = y
    1876         #print_arr("ChoSol: fl", fl)
    1877    
    1878         # next, compute bl from fl and b
    1879         bl[0] = b[0] / fl[0][0]
    1880         for i in (1, 2):
    1881             z = 0.0
    1882             for k in range(i):
    1883                 z += fl[i][k] * bl[k]
    1884                 #print "\t", i, k, z
    1885             bl[i] = (b[i] - z) / fl[i][i]
    1886         #print_vec("ChoSol: bl", bl)
    1887    
    1888         # last, compute beta from bl and fl
    1889         beta = np.empty((n))
    1890         beta[-1] = bl[-1] / fl[-1][-1]
    1891         for i in (1, 0):
    1892             z = 0.0
    1893             for k in range(i+1, n):
    1894                 z += fl[k][i] * beta[k]
    1895                 #print "\t\t", i, k, 'z=', z
    1896             beta[i] = (bl[i] - z) / fl[i][i]
    1897         #print_vec("ChoSol: beta", beta)
    1898    
    1899         return beta
    1900 
    1901     def MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2):
    1902         '''
    1903         move beta one step closer towards the solution
    1904         '''
    1905         a_lower, a_upper = 0., 1.          # bracket  "a"
    1906         cmin, beta = ChiNow (a_lower, c1, c2, s1, s2)
    1907         #print "MaxEntMove: cmin = %g" % cmin
    1908         if cmin*chisq > chizer:
    1909             ctarg = (1.0 + cmin)/2
    1910         else:
    1911             ctarg = chizer/chisq
    1912         f_lower = cmin - ctarg
    1913         c_upper, beta = ChiNow (a_upper, c1, c2, s1, s2)
    1914         f_upper = c_upper - ctarg
    1915    
    1916         fx = 2*MOVE_PASSES      # just to start off
    1917         loop = 1
    1918         while abs(fx) >= MOVE_PASSES and loop <= MAX_MOVE_LOOPS:
    1919             a_new = (a_lower + a_upper) * 0.5           # search by bisection
    1920             c_new, beta = ChiNow (a_new, c1, c2, s1, s2)
    1921             fx = c_new - ctarg
    1922             # tighten the search range for the next pass
    1923             if f_lower*fx > 0:
    1924                 a_lower, f_lower = a_new, fx
    1925             if f_upper*fx > 0:
    1926                 a_upper, f_upper = a_new, fx
    1927             loop += 1
    1928    
    1929         if abs(fx) >= MOVE_PASSES or loop > MAX_MOVE_LOOPS:
    1930             msg = "MaxEntMove: Loop counter = "
    1931             msg += str(MAX_MOVE_LOOPS)
    1932             msg += '  No convergence in alpha chop'
    1933             raise MaxEntException(msg)
    1934    
    1935         w = Dist (s2, beta);
    1936         m = SEARCH_DIRECTIONS
    1937         if (w > DISTANCE_LIMIT_FACTOR*fSum/blank):        # invoke the distance penalty, SB eq. 17
    1938             for k in range(m):
    1939                 beta[k] *= math.sqrt (fSum/(blank*w))
    1940         chtarg = ctarg * chisq
    1941         return w, chtarg, loop, a_new, fx, beta
    1942        
    1943 #MaxEnt_SB starts here   
    1944 
    1945     if image_to_data == None:
    1946         image_to_data = opus
    1947     if data_to_image == None:
    1948         data_to_image = tropus
    1949     n   = len(base)
    1950     npt = len(datum)
    1951 
    1952     # Note that the order of subscripts for
    1953     # "xi" and "eta" has been reversed from
    1954     # the convention used in the FORTRAN version
    1955     # to enable parts of them to be passed as
    1956     # as vectors to "image_to_data" and "data_to_image".
    1957     xi      = np.zeros((SEARCH_DIRECTIONS, n))
    1958     eta     = np.zeros((SEARCH_DIRECTIONS, npt))
    1959     beta    = np.zeros((SEARCH_DIRECTIONS))
    1960     # s1      = np.zeros((SEARCH_DIRECTIONS))
    1961     # c1      = np.zeros((SEARCH_DIRECTIONS))
    1962     s2      = np.zeros((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
    1963     c2      = np.zeros((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
    1964 
    1965     # TODO: replace blank (scalar) with base (vector)
    1966     blank = sum(base) / len(base)   # use the average value of base
    1967 
    1968     chizer, chtarg = npt*1.0, npt*1.0
    1969     f = base * 1.0                              # starting distribution is base
    1970 
    1971     fSum  = sum(f)                              # find the sum of the f-vector
    1972     z = (datum - image_to_data (f, G)) / sigma  # standardized residuals, SB eq. 3
    1973     chisq = sum(z*z)                            # Chi^2, SB eq. 4
    1974 
    1975     for iter in range(IterMax):
    1976         ox = -2 * z / sigma                        # gradient of Chi^2
    1977 
    1978         cgrad = data_to_image (ox, G)              # cgrad[i] = del(C)/del(f[i]), SB eq. 8
    1979         sgrad = -np.log(f/base) / (blank*math.exp (1.0))  # sgrad[i] = del(S)/del(f[i])
    1980         snorm = math.sqrt(sum(f * sgrad*sgrad))    # entropy term, SB eq. 22
    1981         cnorm = math.sqrt(sum(f * cgrad*cgrad))    # ChiSqr term, SB eq. 22
    1982         tnorm = sum(f * sgrad * cgrad)             # norm for gradient term TEST
    1983 
    1984         a = 1.0
    1985         b = 1.0 / cnorm
    1986         if iter == 0:
    1987             test = 0.0     # mismatch between entropy and ChiSquared gradients
    1988         else:
    1989             test = math.sqrt( ( 1.0 - tnorm/(snorm*cnorm) )/2 ) # SB eq. 37?
    1990             a = 0.5 / (snorm * test)
    1991             b *= 0.5 / test
    1992         xi[0] = f * cgrad / cnorm
    1993         xi[1] = f * (a * sgrad - b * cgrad)
    1994 
    1995         eta[0] = image_to_data (xi[0], G);          # image --> data
    1996         eta[1] = image_to_data (xi[1], G);          # image --> data
    1997         ox = eta[1] / (sigma * sigma)
    1998         xi[2] = data_to_image (ox, G);              # data --> image
    1999         a = 1.0 / math.sqrt(sum(f * xi[2]*xi[2]))
    2000         xi[2] = f * xi[2] * a
    2001         eta[2] = image_to_data (xi[2], G)           # image --> data
    2002        
    2003 #         print_arr("MaxEnt: eta.transpose()", eta.transpose())
    2004 #         print_arr("MaxEnt: xi.transpose()", xi.transpose())
    2005 
    2006         # prepare the search directions for the conjugate gradient technique
    2007         c1 = xi.dot(cgrad) / chisq                          # C_mu, SB eq. 24
    2008         s1 = xi.dot(sgrad)                          # S_mu, SB eq. 24
    2009 #         print_vec("MaxEnt: c1", c1)
    2010 #         print_vec("MaxEnt: s1", s1)
    2011 
    2012         for k in range(SEARCH_DIRECTIONS):
    2013             for l in range(k+1):
    2014                 c2[k][l] = 2 * sum(eta[k] * eta[l] / sigma/sigma) / chisq
    2015                 s2[k][l] = -sum(xi[k] * xi[l] / f) / blank
    2016 #         print_arr("MaxEnt: c2", c2)
    2017 #         print_arr("MaxEnt: s2", s2)
    2018 
    2019         # reflect across the body diagonal
    2020         for k, l in ((0,1), (0,2), (1,2)):
    2021             c2[k][l] = c2[l][k]                     #  M_(mu,nu)
    2022             s2[k][l] = s2[l][k]                     #  g_(mu,nu)
    2023  
    2024         beta[0] = -0.5 * c1[0] / c2[0][0]
    2025         beta[1] = 0.0
    2026         beta[2] = 0.0
    2027         if (iter > 0):
    2028             w, chtarg, loop, a_new, fx, beta = MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2)
    2029 
    2030         f_old = f.copy()    # preserve the last image
    2031         f += xi.transpose().dot(beta)   # move the image towards the solution, SB eq. 25
    2032        
    2033         # As mentioned at the top of p.119,
    2034         # need to protect against stray negative values.
    2035         # In this case, set them to RESET_STRAYS * base[i]
    2036         #f = f.clip(RESET_STRAYS * blank, f.max())
    2037         for i in range(n):
    2038             if f[i] <= 0.0:
    2039                 f[i] = RESET_STRAYS * base[i]
    2040         df = f - f_old
    2041         fSum = sum(f)
    2042         fChange = sum(df)
    2043 
    2044         # calculate the normalized entropy
    2045         S = sum((f/fSum) * np.log(f/fSum))      # normalized entropy, S&B eq. 1
    2046         z = (datum - image_to_data (f, G)) / sigma  # standardized residuals
    2047         chisq = sum(z*z)                            # report this ChiSq
    2048 
    2049         if report:
    2050             print " MaxEnt trial/max: %3d/%3d" % ((iter+1), IterMax)
    2051             print " Residual: %5.2lf%% Entropy: %8lg" % (100*test, S)
    2052             if iter > 0:
    2053                 value = 100*( math.sqrt(chisq/chtarg)-1)
    2054             else:
    2055                 value = 0
    2056       #      print " %12.5lg %10.4lf" % ( math.sqrt(chtarg/npt), value )
    2057             print " Function sum: %.6lg Change from last: %.2lf%%\n" % (fSum,100*fChange/fSum)
    2058 
    2059         # See if we have finished our task.
    2060         # do the hardest test first
    2061         if (abs(chisq/chizer-1.0) < CHI_SQR_LIMIT) and  (test < TEST_LIMIT):
    2062             print ' Convergence achieved.'
    2063             return chisq,f,image_to_data(f, G)     # solution FOUND returns here
    2064     print ' No convergence! Try increasing Error multiplier.'
    2065     return chisq,f,image_to_data(f, G)       # no solution after IterMax iterations
    2066 
    2067    
    2068 ###############################################################################
    2069 #### IPG/TNNLS Routines
    2070 ###############################################################################
    2071 
    2072 def IPG(datum,sigma,G,Bins,Dbins,IterMax,Qvec=[],approach=0.8,Power=-1,report=False):
    2073     ''' An implementation of the Interior-Point Gradient method of
    2074     Michael Merritt & Yin Zhang, Technical Report TR04-08, Dept. of Comp. and
    2075     Appl. Math., Rice Univ., Houston, Texas 77005, U.S.A. found on the web at
    2076     http://www.caam.rice.edu/caam/trs/2004/TR04-08.pdf
    2077     Problem addressed: Total Non-Negative Least Squares (TNNLS)
    2078     :param float datum[]:
    2079     :param float sigma[]:
    2080     :param float[][] G: transformation matrix
    2081     :param int IterMax:
    2082     :param float Qvec: data positions for Power = 0-4
    2083     :param float approach: 0.8 default fitting parameter
    2084     :param int Power: 0-4 for Q^Power weighting, -1 to use input sigma
    2085    
    2086     '''
    2087     if Power < 0:
    2088         GmatE = G/sigma[:np.newaxis]
    2089         IntE = datum/sigma
    2090         pwr = 0
    2091         QvecP = np.ones_like(datum)
    2092     else:
    2093         GmatE = G[:]
    2094         IntE = datum[:]
    2095         pwr = Power
    2096         QvecP = Qvec**pwr
    2097     Amat = GmatE*QvecP[:np.newaxis]
    2098     AAmat = np.inner(Amat,Amat)
    2099     Bvec = IntE*QvecP
    2100     Xw = np.ones_like(Bins)*1.e-32
    2101     calc = np.dot(G.T,Xw)
    2102     nIter = 0
    2103     err = 10.
    2104     while (nIter<IterMax) and (err > 1.):
    2105         #Step 1 in M&Z paper:
    2106         Qk = np.inner(AAmat,Xw)-np.inner(Amat,Bvec)
    2107         Dk = Xw/np.inner(AAmat,Xw)
    2108         Pk = -Dk*Qk
    2109         #Step 2 in M&Z paper:
    2110         alpSt = -np.inner(Pk,Qk)/np.inner(Pk,np.inner(AAmat,Pk))
    2111         alpWv = -Xw/Pk
    2112         AkSt = [np.where(Pk[i]<0,np.min((approach*alpWv[i],alpSt)),alpSt) for i in range(Pk.shape[0])]
    2113         #Step 3 in M&Z paper:
    2114         shift = AkSt*Pk
    2115         Xw += shift
    2116         #done IPG; now results
    2117         nIter += 1
    2118         calc = np.dot(G.T,Xw)
    2119         chisq = np.sum(((datum-calc)/sigma)**2)
    2120         err = chisq/len(datum)
    2121         if report:
    2122             print ' Iteration: %d, chisq: %.3g, sum(shift^2): %.3g'%(nIter,chisq,np.sum(shift**2))
    2123     return chisq,Xw,calc
    2124 
    2125 ###############################################################################
    2126 #### SASD Utilities
    2127 ###############################################################################
    2128 
    2129 def SetScale(Data,refData):
    2130     Profile,Limits,Sample = Data
    2131     refProfile,refLimits,refSample = refData
    2132     x,y = Profile[:2]
    2133     rx,ry = refProfile[:2]
    2134     Beg = np.max([rx[0],x[0],Limits[1][0],refLimits[1][0]])
    2135     Fin = np.min([rx[-1],x[-1],Limits[1][1],refLimits[1][1]])
    2136     iBeg = np.searchsorted(x,Beg)
    2137     iFin = np.searchsorted(x,Fin)
    2138     sum = np.sum(y[iBeg:iFin])
    2139     refsum = np.sum(np.interp(x[iBeg:iFin],rx,ry,0,0))
    2140     Sample['Scale'][0] = refSample['Scale'][0]*refsum/sum
    2141    
    2142 def Bestimate(G,Rg,P):
    2143     return (G*P/Rg**P)*np.exp(scsp.gammaln(P/2))
    2144    
    2145 ###############################################################################
    2146 #### Size distribution
    2147 ###############################################################################
    2148 
    2149 def SizeDistribution(Profile,ProfDict,Limits,Substances,Sample,data):
    2150     shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
    2151         'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
    2152         'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
    2153         'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
    2154         'Cylinder diam':[CylinderDFF,CylinderDVol]}
    2155     Shape = data['Size']['Shape'][0]
    2156     Parms = data['Size']['Shape'][1:]
    2157     if data['Size']['logBins']:
    2158         Bins = np.logspace(np.log10(data['Size']['MinDiam']),np.log10(data['Size']['MaxDiam']),
    2159             data['Size']['Nbins']+1,True)/2.        #make radii
    2160     else:
    2161         Bins = np.linspace(data['Size']['MinDiam'],data['Size']['MaxDiam'],
    2162             data['Size']['Nbins']+1,True)/2.        #make radii
    2163     Dbins = np.diff(Bins)
    2164     Bins = Bins[:-1]+Dbins/2.
    2165     Contrast = Sample['Contrast'][1]
    2166     Scale = Sample['Scale'][0]
    2167     Sky = 10**data['Size']['MaxEnt']['Sky']
    2168     BinsBack = np.ones_like(Bins)*Sky*Scale/Contrast
    2169     Back = data['Back']
    2170     Q,Io,wt,Ic,Ib = Profile[:5]
    2171     Qmin = Limits[1][0]
    2172     Qmax = Limits[1][1]
    2173     wtFactor = ProfDict['wtFactor']
    2174     Ibeg = np.searchsorted(Q,Qmin)
    2175     Ifin = np.searchsorted(Q,Qmax)
    2176     BinMag = np.zeros_like(Bins)
    2177     Ic[:] = 0.
    2178     Gmat = G_matrix(Q[Ibeg:Ifin],Bins,Contrast,shapes[Shape][0],shapes[Shape][1],args=Parms)
    2179     if 'MaxEnt' == data['Size']['Method']:
    2180         chisq,BinMag,Ic[Ibeg:Ifin] = MaxEnt_SB(Scale*Io[Ibeg:Ifin]-Back[0],
    2181             Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),Gmat,BinsBack,
    2182             data['Size']['MaxEnt']['Niter'],report=True)
    2183     elif 'IPG' == data['Size']['Method']:
    2184         chisq,BinMag,Ic[Ibeg:Ifin] = IPG(Scale*Io[Ibeg:Ifin]-Back[0],Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),
    2185             Gmat,Bins,Dbins,data['Size']['IPG']['Niter'],Q[Ibeg:Ifin],approach=0.8,
    2186             Power=data['Size']['IPG']['Power'],report=True)
    2187     Ib[:] = Back[0]
    2188     Ic[Ibeg:Ifin] += Back[0]
    2189     print ' Final chi^2: %.3f'%(chisq)
    2190     Vols = shapes[Shape][1](Bins,Parms)
    2191     data['Size']['Distribution'] = [Bins,Dbins,BinMag/(2.*Dbins)]
    2192        
    2193 ################################################################################
    2194 #### Modelling
    2195 ################################################################################
    2196 
    2197 def ModelFit(Profile,ProfDict,Limits,Substances,Sample,Model):
    2198     shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
    2199         'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
    2200         'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
    2201         'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
    2202         'Unified tube':[UniTubeFF,UniTubeVol],'Cylinder diam':[CylinderDFF,CylinderDVol]}
    2203 
    2204     def GetModelParms():
    2205         parmOrder = ['Volume','Radius','Mean','StdDev','MinSize','G','Rg','B','P','Cutoff',
    2206             'PkInt','PkPos','PkSig','PkGam',]
    2207         parmDict = {}
    2208         varyList = []
    2209         values = []
    2210         levelTypes = []
    2211         Back = Model['Back']
    2212         if Back[1]:
    2213             varyList += ['Back',]
    2214             values.append(Back[0])
    2215         parmDict['Back'] = Back[0]
    2216         partData = Model['Particle']
    2217         parmDict['Matrix density'] = Substances['Substances'][partData['Matrix']['Name']].get('XAnom density',0.0)
    2218         for i,level in enumerate(partData['Levels']):
    2219             cid = str(i)+':'
    2220             controls = level['Controls']
    2221             Type = controls['DistType']
    2222             levelTypes.append(Type)
    2223             if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
    2224                 if Type not in ['Monodosperse',]:
    2225                     parmDict[cid+'NumPoints'] = controls['NumPoints']
    2226                     parmDict[cid+'Cutoff'] = controls['Cutoff']
    2227                 parmDict[cid+'FormFact'] = shapes[controls['FormFact']][0]
    2228                 parmDict[cid+'FFVolume'] = shapes[controls['FormFact']][1]
    2229                 parmDict[cid+'XAnom density'] = Substances['Substances'][controls['Material']].get('XAnom density',0.0)
    2230                 for item in ['Aspect ratio','Length','Thickness','Diameter',]:
    2231                     if item in controls['FFargs']:
    2232                         parmDict[cid+item] = controls['FFargs'][item][0]
    2233                         if controls['FFargs'][item][1]:
    2234                             varyList.append(cid+item)
    2235                             values.append(controls['FFargs'][item][0])
    2236             distDict = controls['DistType']
    2237             for item in parmOrder:
    2238                 if item in level[distDict]:
    2239                     parmDict[cid+item] = level[distDict][item][0]
    2240                     if level[distDict][item][1]:
    2241                         values.append(level[distDict][item][0])
    2242                         varyList.append(cid+item)
    2243         return levelTypes,parmDict,varyList,values
    2244        
    2245     def SetModelParms():
    2246         print ' Refined parameters:'
    2247         if 'Back' in varyList:
    2248             Model['Back'][0] = parmDict['Back']
    2249             print '  %15s %15.4f esd: %15.4g'%('Background:',parmDict['Back'],sigDict['Back'])
    2250         partData = Model['Particle']
    2251         for i,level in enumerate(partData['Levels']):
    2252             Type = level['Controls']['DistType']
    2253             print ' Component %d: Type: %s:'%(i,Type)
    2254             cid = str(i)+':'
    2255             controls = level['Controls']
    2256             Type = controls['DistType']
    2257             if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
    2258                 for item in ['Aspect ratio','Length','Thickness','Diameter',]:
    2259                     if cid+item in varyList:
    2260                         controls['FFargs'][item][0] = parmDict[cid+item]
    2261                         print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
    2262             distDict = controls['DistType']
    2263             for item in level[distDict]:
    2264                 if cid+item in varyList:
    2265                     level[distDict][item][0] = parmDict[cid+item]
    2266                     print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
    2267                    
    2268     def calcSASD(values,Q,Io,wt,levelTypes,parmDict,varyList):
    2269         parmDict.update(zip(varyList,values))
    2270         M = np.sqrt(wt)*(getSASD(Q,levelTypes,parmDict)-Io)
    2271         return M
    2272        
    2273     def getSASD(Q,levelTypes,parmDict):
    2274         Ic = np.ones_like(Q)
    2275         Ic *= parmDict['Back']
    2276         rhoMat = parmDict['Matrix density']
    2277         for i,Type in enumerate(levelTypes):
    2278             cid = str(i)+':'
    2279             if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm']:
    2280                 FFfxn = parmDict[cid+'FormFact']
    2281                 Volfxn = parmDict[cid+'FFVolume']
    2282                 FFargs = []
    2283                 for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter',]:
    2284                     if item in parmDict:
    2285                         FFargs.append(parmDict[item])
    2286                 distDict = {}
    2287                 for item in [cid+'Volume',cid+'Mean',cid+'StdDev',cid+'MinSize',]:
    2288                     if item in parmDict:
    2289                         distDict[item.split(':')[1]] = parmDict[item]
    2290                 rho = parmDict[cid+'XAnom density']
    2291                 contrast = rho**2-rhoMat**2
    2292                 rBins,dBins,dist = MakeDiamDist(Type,parmDict[cid+'NumPoints'],parmDict[cid+'Cutoff'],distDict)
    2293                 Gmat = G_matrix(Q,rBins,contrast,FFfxn,Volfxn,FFargs).T
    2294                 dist *= parmDict[cid+'Volume']
    2295                 Ic += np.dot(Gmat,dist)
    2296             elif 'Unified' in Type:
    2297                 Rg,G,B,P,Rgco = parmDict[cid+'Rg'],parmDict[cid+'G'],parmDict[cid+'B'], \
    2298                     parmDict[cid+'P'],parmDict[cid+'Cutoff']
    2299                 Qstar = Q/(scsp.erf(Q*Rg/np.sqrt(6)))**3
    2300                 Guin = G*np.exp(-(Q*Rg)**2/3)
    2301                 Porod = (B/Qstar**P)*np.exp(-(Q*Rgco)**2/3)
    2302                 Ic += Guin+Porod
    2303             elif 'Porod' in Type:
    2304                 B,P,Rgco = parmDict[cid+'B'],parmDict[cid+'P'],parmDict[cid+'Cutoff']
    2305                 Porod = (B/Q**P)*np.exp(-(Q*Rgco)**2/3)
    2306                 Ic += Porod
    2307             elif 'Mono' in Type:
    2308                 FFfxn = parmDict[cid+'FormFact']
    2309                 Volfxn = parmDict[cid+'FFVolume']
    2310                 FFargs = []
    2311                 for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter',]:
    2312                     if item in parmDict:
    2313                         FFargs.append(parmDict[item])
    2314                 rho = parmDict[cid+'XAnom density']
    2315                 contrast = rho**2-rhoMat**2
    2316                 R = parmDict[cid+'Radius']
    2317                 Gmat = G_matrix(Q,R,contrast,FFfxn,Volfxn,FFargs)             
    2318                 Ic += Gmat[0]*parmDict[cid+'Volume']
    2319             elif 'Bragg' in distFxn:
    2320                 Ic += parmDict[cid+'PkInt']*G2pwd.getPsVoigt(parmDict[cid+'PkPos'],
    2321                     parmDict[cid+'PkSig'],parmDict[cid+'PkGam'],Q)
    2322         return Ic
    2323        
    2324     Q,Io,wt,Ic,Ib = Profile[:5]
    2325     Qmin = Limits[1][0]
    2326     Qmax = Limits[1][1]
    2327     wtFactor = ProfDict['wtFactor']
    2328     Ibeg = np.searchsorted(Q,Qmin)
    2329     Ifin = np.searchsorted(Q,Qmax)
    2330     Ic[:] = 0
    2331     levelTypes,parmDict,varyList,values = GetModelParms()
    2332     result = so.leastsq(calcSASD,values,full_output=True,   #ftol=Ftol,
    2333         args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],levelTypes,parmDict,varyList))
    2334     ncyc = int(result[2]['nfev']/2)
    2335     chisq = np.sum(result[2]['fvec']**2)
    2336     Rvals = {}
    2337     Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
    2338     Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
    2339     parmDict.update(zip(varyList,result[0]))
    2340     Ic[Ibeg:Ifin] = getSASD(Q[Ibeg:Ifin],levelTypes,parmDict)
    2341     try:
    2342         sig = np.sqrt(np.diag(result[1])*Rvals['GOF'])
    2343         sigDict = dict(zip(varyList,sig))
    2344         print ' Results of small angle data modelling fit:'
    2345         print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',Ifin-Ibeg,' Number of parameters: ',len(varyList)
    2346         print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
    2347         SetModelParms()
    2348         covMatrix = result[1]*Rvals['GOF']
    2349         return True,result,varyList,sig,Rvals,covMatrix
    2350     except ValueError:
    2351         return False,0,0,0,0,0
    2352    
    2353 def ModelFxn(Profile,ProfDict,Limits,Substances,Sample,sasdData):
    2354     shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
    2355         'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
    2356         'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
    2357         'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
    2358         'Unified tube':[UniTubeFF,UniTubeVol],'Cylinder diam':[CylinderDFF,CylinderDVol]}
    2359 #    pdb.set_trace()
    2360     partData = sasdData['Particle']
    2361     rhoMat = Substances['Substances'][partData['Matrix']['Name']].get('XAnom density',0.0)
    2362     matFrac = partData['Matrix']['VolFrac'] #[value,flag]       
    2363     Scale = Sample['Scale']     #[value,flag]
    2364     Back = sasdData['Back']
    2365     Q,Io,wt,Ic,Ib = Profile[:5]
    2366     Qmin = Limits[1][0]
    2367     Qmax = Limits[1][1]
    2368     wtFactor = ProfDict['wtFactor']
    2369     Ibeg = np.searchsorted(Q,Qmin)
    2370     Ifin = np.searchsorted(Q,Qmax)
    2371     Ib[:] = Back[0]
    2372     Ic[:] = 0
    2373     Ic[Ibeg:Ifin] += Back[0]
    2374     Rbins = []
    2375     Dist = []
    2376     for level in partData['Levels']:
    2377         controls = level['Controls']
    2378         distFxn = controls['DistType']
    2379         if distFxn in ['LogNormal','Gaussian','LSW','Schulz-Zimm']:
    2380             FFfxn = shapes[controls['FormFact']][0]
    2381             Volfxn = shapes[controls['FormFact']][1]
    2382             FFargs = []
    2383             for item in ['Aspect ratio','Length','Thickness','Diameter',]:
    2384                 if item in controls['FFargs']:
    2385                     FFargs.append(controls['FFargs'][item][0])
    2386             rho = Substances['Substances'][level['Controls']['Material']].get('XAnom density',0.0)
    2387             contrast = rho**2-rhoMat**2
    2388             parmDict = level[controls['DistType']]
    2389             distDict = {}
    2390             for item in parmDict:
    2391                 distDict[item] = parmDict[item][0]
    2392             rBins,dBins,dist = MakeDiamDist(controls['DistType'],controls['NumPoints'],controls['Cutoff'],distDict)
    2393             Gmat = G_matrix(Q[Ibeg:Ifin],rBins,contrast,FFfxn,Volfxn,FFargs).T
    2394             dist *= level[distFxn]['Volume'][0]
    2395             Ic[Ibeg:Ifin] += np.dot(Gmat,dist)
    2396             Rbins.append(rBins)
    2397             Dist.append(dist/(4.*dBins))
    2398         elif 'Unified' in distFxn:
    2399             parmDict = level[controls['DistType']]
    2400             Rg,G,B,P,Rgco = parmDict['Rg'][0],parmDict['G'][0],parmDict['B'][0],    \
    2401                 parmDict['P'][0],parmDict['Cutoff'][0]
    2402             Qstar = Q[Ibeg:Ifin]/(scsp.erf(Q[Ibeg:Ifin]*Rg/np.sqrt(6)))**3
    2403             Guin = G*np.exp(-(Q[Ibeg:Ifin]*Rg)**2/3)
    2404             Porod = (B/Qstar**P)*np.exp(-(Q[Ibeg:Ifin]*Rgco)**2/3)
    2405             Ic[Ibeg:Ifin] += Guin+Porod
    2406             Rbins.append([])
    2407             Dist.append([])
    2408         elif 'Porod' in distFxn:
    2409             parmDict = level[controls['DistType']]
    2410             B,P,Rgco = parmDict['B'][0],parmDict['P'][0],parmDict['Cutoff'][0]
    2411             Porod = (B/Q[Ibeg:Ifin]**P)*np.exp(-(Q[Ibeg:Ifin]*Rgco)**2/3)
    2412             Ic[Ibeg:Ifin] += Porod
    2413             Rbins.append([])
    2414             Dist.append([])
    2415         elif 'Mono' in distFxn:
    2416             FFfxn = shapes[controls['FormFact']][0]
    2417             Volfxn = shapes[controls['FormFact']][1]
    2418             FFargs = []
    2419             for item in ['Aspect ratio','Length','Thickness','Diameter',]:
    2420                 if item in controls['FFargs']:
    2421                     FFargs.append(controls['FFargs'][item][0])
    2422             rho = Substances['Substances'][level['Controls']['Material']].get('XAnom density',0.0)
    2423             contrast = rho**2-rhoMat**2
    2424             R = level[controls['DistType']]['Radius'][0]
    2425             Gmat = G_matrix(Q[Ibeg:Ifin],R,contrast,FFfxn,Volfxn,FFargs)             
    2426             Ic[Ibeg:Ifin] += Gmat[0]*level[distFxn]['Volume'][0]
    2427             Rbins.append([])
    2428             Dist.append([])
    2429         elif 'Bragg' in distFxn:
    2430             parmDict = level[controls['DistType']]
    2431             Ic[Ibeg:Ifin] += parmDict['PkInt'][0]*G2pwd.getPsVoigt(parmDict['PkPos'][0],
    2432                 parmDict['PkSig'][0],parmDict['PkGam'][0],Q[Ibeg:Ifin])
    2433             Rbins.append([])
    2434             Dist.append([])           
    2435     sasdData['Size Calc'] = [Rbins,Dist]
    2436    
    2437 def MakeDiamDist(DistName,nPoints,cutoff,distDict):
    2438    
    2439     if 'LogNormal' in DistName:
    2440         distFxn = 'LogNormalDist'
    2441         cumeFxn = 'LogNormalCume'
    2442         pos = distDict['MinSize']
    2443         args = [distDict['Mean'],distDict['StdDev']]
    2444         step = 4.*np.sqrt(np.exp(distDict['StdDev']**2)*(np.exp(distDict['StdDev']**2)-1.))
    2445         mode = distDict['MinSize']+distDict['Mean']/np.exp(distDict['StdDev']**2)
    2446         minX = 1. #pos
    2447         maxX = 1.e4 #np.min([mode+nPoints*step*40,1.e5])
    2448     elif 'Gauss' in DistName:
    2449         distFxn = 'GaussDist'
    2450         cumeFxn = 'GaussCume'
    2451         pos = distDict['Mean']
    2452         args = [distDict['StdDev']]
    2453         step = 0.02*distDict['StdDev']
    2454         mode = distDict['Mean']
    2455         minX = np.max([mode-4.*distDict['StdDev'],1.])
    2456         maxX = np.min([mode+4.*distDict['StdDev'],1.e5])
    2457     elif 'LSW' in DistName:
    2458         distFxn = 'LSWDist'
    2459         cumeFxn = 'LSWCume'
    2460         pos = distDict['Mean']
    2461         args = []
    2462         minX,maxX = [0.,2.*pos]
    2463     elif 'Schulz' in DistName:
    2464         distFxn = 'SchulzZimmDist'
    2465         cumeFxn = 'SchulzZimmCume'
    2466         pos = distDict['Mean']
    2467         args = [distDict['StdDev']]
    2468         minX = np.max([1.,pos-4.*distDict['StdDev']])
    2469         maxX = np.min([pos+4.*distDict['StdDev'],1.e5])
    2470     nP = 500
    2471     Diam = np.logspace(0.,5.,nP,True)
    2472     TCW = eval(cumeFxn+'(Diam,pos,args)')
    2473     CumeTgts = np.linspace(cutoff,(1.-cutoff),nPoints+1,True)
    2474     rBins = np.interp(CumeTgts,TCW,Diam,0,0)
    2475     dBins = np.diff(rBins)
    2476     rBins = rBins[:-1]+dBins/2.
    2477     return rBins,dBins,eval(distFxn+'(rBins,pos,args)')
    2478 
    2479 ################################################################################
    2480 #### MaxEnt testing stuff
    2481 ################################################################################
    2482 
    2483 def print_vec(text, a):
    2484     '''print the contents of a vector to the console'''
    2485     n = a.shape[0]
    2486     print "%s[ = (" % text,
    2487     for i in range(n):
    2488         s = " %g, " % a[i]
    2489         print s,
    2490     print ")"
    2491 
    2492 def print_arr(text, a):
    2493     '''print the contents of an array to the console'''
    2494     n, m = a.shape
    2495     print "%s[][] = (" % text
    2496     for i in range(n):
    2497         print " (",
    2498         for j in range(m):
    2499             print " %g, " % a[i][j],
    2500         print "),"
    2501     print ")"
    2502 
    2503 def test_MaxEnt_SB(report=True):
    2504     def readTextData(filename):
    2505         '''return q, I, dI from a 3-column text file'''
    2506         if not os.path.exists(filename):
    2507             raise Exception("file not found: " + filename)
    2508         buf = [line.split() for line in open(filename, 'r').readlines()]
    2509         M = len(buf)
    2510         buf = zip(*buf)         # transpose rows and columns
    2511         q  = np.array(buf[0], dtype=np.float64)
    2512         I  = np.array(buf[1], dtype=np.float64)
    2513         dI = np.array(buf[2], dtype=np.float64)
    2514         return q, I, dI
    2515     print "MaxEnt_SB: "
    2516     test_data_file = os.path.join( 'testinp', 'test.sas')
    2517     rhosq = 100     # scattering contrast, 10^20 1/cm^-4
    2518     bkg   = 0.1     #   I = I - bkg
    2519     dMin, dMax, nRadii = 25, 9000, 40
    2520     defaultDistLevel = 1.0e-6
    2521     IterMax = 40
    2522     errFac = 1.05
    2523    
    2524     r    = np.logspace(math.log10(dMin), math.log10(dMax), nRadii)/2
    2525     dr   = r * (r[1]/r[0] - 1)          # step size
    2526     f_dr = np.ndarray((nRadii)) * 0  # volume fraction histogram
    2527     b    = np.ndarray((nRadii)) * 0 + defaultDistLevel  # MaxEnt "sky background"
    2528    
    2529     qVec, I, dI = readTextData(test_data_file)
    2530     G = G_matrix(qVec,r,rhosq,SphereFF,SphereVol,args=())
    2531    
    2532     chisq,f_dr,Ic = MaxEnt_SB(I - bkg, dI*errFac, b, IterMax, G, report=report)
    2533     if f_dr is None:
    2534         print "no solution"
    2535         return
    2536    
    2537     print "solution reached"
    2538     for a,b,c in zip(r.tolist(), dr.tolist(), f_dr.tolist()):
    2539         print '%10.4f %10.4f %12.4g'%(a,b,c)
    2540 
    2541 def tests():
    2542     test_MaxEnt_SB(report=True)
    2543 
    2544 if __name__ == '__main__':
    2545     tests()
    2546 
    2547 >>>>>>> .r1299
Note: See TracChangeset for help on using the changeset viewer.