Changeset 1302 for trunk/GSASIIsasd.py
 Timestamp:
 Apr 25, 2014 12:56:35 PM (7 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIsasd.py
r1301 r1302 1 <<<<<<< .mine2 1 #/usr/bin/env python 3 2 # * coding: utf8 * … … 1278 1277 tests() 1279 1278 1280 =======1281 #/usr/bin/env python1282 # * coding: utf8 *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 os1296 import sys1297 import math1298 import time1299 1300 import numpy as np1301 import scipy as sp1302 import numpy.linalg as nl1303 from numpy.fft import ifft, fft, fftshift1304 import scipy.special as scsp1305 import scipy.interpolate as si1306 import scipy.stats as st1307 import scipy.optimize as so1308 #import pdb1309 1310 import GSASIIpath1311 GSASIIpath.SetVersionNumber("$Revision$")1312 import GSASIIlattice as G2lat1313 import GSASIIspc as G2spc1314 import GSASIIElem as G2elem1315 import GSASIIgrid as G2gd1316 import GSASIIIO as G2IO1317 import GSASIImath as G2mth1318 import GSASIIpwd as G2pwd1319 1320 # trig functions in degrees1321 sind = lambda x: math.sin(x*math.pi/180.)1322 asind = lambda x: 180.*math.asin(x)/math.pi1323 tand = lambda x: math.tan(x*math.pi/180.)1324 atand = lambda x: 180.*math.atan(x)/math.pi1325 atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi1326 cosd = lambda x: math.cos(x*math.pi/180.)1327 acosd = lambda x: 180.*math.acos(x)/math.pi1328 rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)1329 #numpy versions1330 npsind = lambda x: np.sin(x*np.pi/180.)1331 npasind = lambda x: 180.*np.arcsin(x)/math.pi1332 npcosd = lambda x: np.cos(x*math.pi/180.)1333 npacosd = lambda x: 180.*np.arccos(x)/math.pi1334 nptand = lambda x: np.tan(x*math.pi/180.)1335 npatand = lambda x: 180.*np.arctan(x)/np.pi1336 npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi1337 npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave1338 npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)1339 1340 ###############################################################################1341 #### Particle form factors1342 ###############################################################################1343 1344 def SphereFF(Q,R,args=()):1345 ''' Compute hard sphere form factor  can use numpy arrays1346 param float Q: Q value array (usually in A1)1347 param float R: sphere radius (Usually in A  must match Q1 units)1348 param array args: ignored1349 returns float: form factors as array as needed1350 '''1351 QR = Q[:,np.newaxis]*R1352 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 array1357 param float Q : Q value array (usually in A1)1358 param float R: radius along 2 axes of spheroid1359 param array args: [float AR]: aspect ratio so 3rd axis = R*AR1360 returns float: form factors as array as needed1361 '''1362 NP = 501363 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**21.)*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 arrays1373 param float Q: Q value array (A1)1374 param float R: cylinder radius (A)1375 param array args: [float L]: cylinder length (A)1376 returns float: form factor1377 '''1378 L = args[0]1379 NP = 2001380 alp = np.linspace(0,np.pi/2.,NP)1381 QL = Q[:,np.newaxis]*L1382 LBessArg = 0.5*QL[:,:,np.newaxis]**np.cos(alp)1383 LBess = np.where(LBessArg<1.e6,1.,np.sin(LBessArg)/LBessArg)1384 QR = Q[:,np.newaxis]*R1385 SBessArg = QR[:,:,np.newaxis]*np.sin(alp)1386 SBess = np.where(SBessArg<1.e6,0.5,scsp.jv(1,SBessArg)/SBessArg)1387 LSBess = LBess*SBess1388 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 arrays1392 param float Q: Q value array (A1)1393 param float L: cylinder half length (A)1394 param array args: [float R]: cylinder radius (A)1395 returns float: form factor1396 '''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 arrays1402 param float Q: Q value array (A1)1403 param float R: cylinder radius (A)1404 param array args: [float AR]: cylinder aspect ratio = L/D = L/2R1405 returns float: form factor1406 '''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 arrays1412 param float Q: Q value array (A1)1413 param float R: cylinder radius (A)1414 param array args: ignored1415 returns float: form factor1416 '''1417 Rg = np.sqrt(3./5.)*R1418 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)))**31420 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 arrays1424 param float Q: Q value array (A1)1425 param float R: cylinder radius (A)1426 param array args: [float R]: cylinder radius (A)1427 returns float: form factor1428 '''1429 L = args[0]1430 Rg2 = np.sqrt(R**2/2+L**2/12)1431 B2 = np.pi/L1432 Rg1 = np.sqrt(3.)*R/2.1433 G1 = (2./3.)*R/L1434 B1 = 4.*(L+R)/(R**3*L**2)1435 QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg2/np.sqrt(6)))**31436 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)))**31438 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 arrays1443 param float Q: Q value array (A1)1444 param float R: cylinder radius (A)1445 param array args: [float AR]: cylinder aspect ratio = L/D = L/2R1446 returns float: form factor1447 '''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 arrays1453 param float Q: Q value array (A1)1454 param float R: cylinder radius (A)1455 param array args: [float T]: disk thickness (A)1456 returns float: form factor1457 '''1458 T = args[0]1459 Rg2 = np.sqrt(R**2/2.+T**2/12.)1460 B2 = 2./R**21461 Rg1 = np.sqrt(3.)*T/2.1462 RgC2 = 1.1*Rg11463 G1 = (2./3.)*(T/R)**21464 B1 = 4.*(T+R)/(R**3*T**2)1465 QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg2/np.sqrt(6)))**31466 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)))**31468 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 arrays1473 assumes that core of tube is same as the matrix/solvent so contrast1474 is from tube wall vs matrix1475 param float Q: Q value array (A1)1476 param float R: cylinder radius (A)1477 param array args: [float L,T]: tube length & wall thickness(A)1478 returns float: form factor1479 '''1480 L,T = args[:2]1481 Ri = RT1482 DR2 = R**2Ri**21483 Vt = np.pi*DR2*L1484 Rg3 = np.sqrt(DR2/2.+L**2/12.)1485 B1 = 4.*np.pi**2*(DR2+L*(R+Ri))/Vt**21486 B2 = np.pi**2*T/Vt1487 B3 = np.pi/L1488 QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg3/np.sqrt(6)))**31489 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)))**31491 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)))**31493 FF += B1/QstV**41494 return np.sqrt(FF)1495 1496 ###############################################################################1497 #### Particle volumes1498 ###############################################################################1499 1500 def SphereVol(R,args=()):1501 ''' Compute volume of sphere1502  numpy array friendly1503 param float R: sphere radius1504 param array args: ignored1505 returns float: volume1506 '''1507 return (4./3.)*np.pi*R**31508 1509 def SpheroidVol(R,args):1510 ''' Compute volume of cylindrically symmetric ellipsoid (spheroid)1511  numpy array friendly1512 param float R: radius along 2 axes of spheroid1513 param array args: [float AR]: aspect ratio so radius of 3rd axis = R*AR1514 returns float: volume1515 '''1516 AR = args[0]1517 return AR*SphereVol(R)1518 1519 def CylinderVol(R,args):1520 ''' Compute cylinder volume for radius & length1521  numpy array friendly1522 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**21528 1529 def CylinderDVol(L,args):1530 ''' Compute cylinder volume for length & diameter1531  numpy array friendly1532 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/D1541  numpy array friendly1542 param float: R radius (A)1543 param array args: [float AR]: =L/D=L/2R aspect ratio1544 returns float:volume1545 '''1546 AR = args[0]1547 return CylinderVol(R,[2.*R*AR,])1548 1549 def UniSphereVol(R,args=()):1550 ''' Compute volume of sphere1551  numpy array friendly1552 param float R: sphere radius1553 param array args: ignored1554 returns float: volume1555 '''1556 return SphereVol(R)1557 1558 def UniRodVol(R,args):1559 ''' Compute cylinder volume for radius & length1560  numpy array friendly1561 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 ratio1570  numpy array friendly1571 param float R: diameter (A)1572 param array args: [float AR]: =L/D=L/2R aspect ratio1573 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 & thickness1580  numpy array friendly1581 param float R: diameter (A)1582 param array args: [float T]: thickness1583 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 thickness1590  numpy array friendly1591 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 wall1594 '''1595 L,T = args[:2]1596 return CylinderVol(R,[L,])CylinderVol(RT,[L,])1597 1598 ################################################################################1599 #### Distribution functions & their cumulative fxns1600 ################################################################################1601 1602 def LogNormalDist(x,pos,args):1603 ''' Standard LogNormal distribution  numpy friendly on x axis1604 ref: http://www.itl.nist.gov/div898/handbook/index.htm 1.3.6.6.91605 param float x: independent axis (can be numpy array)1606 param float pos: location of distribution1607 param float scale: width of distribution (m)1608 param float shape: shape  (sigma of log(LogNormal))1609 returns float: LogNormal distribution1610 '''1611 scale,shape = args1612 return np.exp(np.log((xpos)/scale)**2/(2.*shape**2))/(np.sqrt(2.*np.pi)*(xpos)*shape)1613 1614 def GaussDist(x,pos,args):1615 ''' Standard Normal distribution  numpy friendly on x axis1616 param float x: independent axis (can be numpy array)1617 param float pos: location of distribution1618 param float scale: width of distribution (sigma)1619 param float shape: not used1620 returns float: Normal distribution1621 '''1622 scale = args[0]1623 return (1./(scale*np.sqrt(2.*np.pi)))*np.exp((xpos)**2/(2.*scale**2))1624 1625 def LSWDist(x,pos,args=[]):1626 ''' LifshitzSlyozovWagner Ostwald ripening distribution  numpy friendly on x axis1627 ref:1628 param float x: independent axis (can be numpy array)1629 param float pos: location of distribution1630 param float scale: not used1631 param float shape: not used1632 returns float: LSW distribution1633 '''1634 redX = x/pos1635 result = (81.*2**(5/3.))*(redX**2*np.exp(redX/(1.5redX)))/((1.5redX)**(11/3.)*(3.redX)**(7/3.))1636 1637 return np.nan_to_num(result/pos)1638 1639 def SchulzZimmDist(x,pos,args):1640 ''' SchulzZimm macromolecule distribution  numpy friendly on x axis1641 ref: http://goldbook.iupac.org/S05502.html1642 param float x: independent axis (can be numpy array)1643 param float pos: location of distribution1644 param float scale: width of distribution (sigma)1645 param float shape: not used1646 returns float: SchulzZimm distribution1647 '''1648 scale = args[0]1649 b = (2.*pos/scale)**21650 a = b/pos1651 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 axis1658 ref: http://www.itl.nist.gov/div898/handbook/index.htm 1.3.6.6.91659 param float x: independent axis (can be numpy array)1660 param float pos: location of distribution1661 param float scale: width of distribution (sigma)1662 param float shape: shape parameter1663 returns float: LogNormal cumulative distribution1664 '''1665 scale,shape = args1666 lredX = np.log((xpos)/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 axis1671 param float x: independent axis (can be numpy array)1672 param float pos: location of distribution1673 param float scale: width of distribution (sigma)1674 param float shape: not used1675 returns float: Normal cumulative distribution1676 '''1677 scale = args[0]1678 redX = (xpos)/scale1679 return (scsp.erf(redX/np.sqrt(2.))+1.)/2.1680 1681 def LSWCume(x,pos,args=[]):1682 ''' LifshitzSlyozovWagner Ostwald ripening cumulative distribution  numpy friendly on x axis1683 param float x: independent axis (can be numpy array)1684 param float pos: location of distribution1685 param float scale: not used1686 param float shape: not used1687 returns float: LSW cumulative distribution1688 '''1689 nP = 5001690 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 ''' SchulzZimm cumulative distribution  numpy friendly on x axis1699 param float x: independent axis (can be numpy array)1700 param float pos: location of distribution1701 param float scale: width of distribution (sigma)1702 param float shape: not used1703 returns float: Normal distribution1704 '''1705 scale = args[0]1706 nP = 5001707 xMin = np.max([0.,pos4.*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 ##### SBMaxEnt1720 ################################################################################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 contrast1728 :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.e4*(contrast*Vol*FF**2).T #10^20 vs 10^241735 1736 '''1737 sbmaxent1738 1739 Entropy maximization routine as described in the article1740 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, LLC1744 :license: This file is distributed subject to a Software License Agreement found1745 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. Workshop1751 Neutron Scattering Data Analysis, Rutherford1752 Appleton Laboratory, UK, 1986; ed. MW Johnson,1753 IOP Conference Series 81 (1986) 81  86, Institute1754 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 pass1766 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 Bryan1770 1771 :param float datum[]:1772 :param float sigma[]:1773 :param float[][] G: transformation matrix1774 :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 convergence1783 CHI_SQR_LIMIT = 0.01 # maximum difference in ChiSqr for a solution1784 SEARCH_DIRECTIONS = 3 # <10. This code requires value = 31785 RESET_STRAYS = 1 # was 0.001, correction of stray negative values1786 DISTANCE_LIMIT_FACTOR = 0.1 # limitation on df to constrain runaways1787 1788 MAX_MOVE_LOOPS = 500 # for no solution in routine: move,1789 MOVE_PASSES = 0.001 # convergence test in routine: move1790 1791 def tropus (data, G):1792 '''1793 tropus: transform dataspace > solutionspace: [G] * data1794 1795 default definition, caller can use this definition or provide an alternative1796 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 solutionspace > dataspace: [G]^tr * image1806 1807 default definition, caller can use this definition or provide an alternative1808 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 = 01818 n = beta.shape[0]1819 for k in range(n):1820 z = sum(s2[k] * beta)1821 w += beta[k] * z1822 return w1823 1824 def ChiNow(ax, c1, c2, s1, s2):1825 '''1826 ChiNow1827 1828 :returns tuple: (ChiNow computation of ``w``, beta)1829 '''1830 1831 bx = 1  ax1832 a = bx * c2  ax * s21833 b = (bx * c1  ax * s1)1834 1835 beta = ChoSol(a, b)1836 w = 1.01837 for k in range(SEARCH_DIRECTIONS):1838 w += beta[k] * (c1[k] + 0.5*sum(c2[k] * beta))1839 return w, beta1840 1841 def ChoSol(a, b):1842 '''1843 ChoSol: ? chop the solution vectors ?1844 1845 :returns: new vector beta1846 '''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 a1861 # note fl is a lower triangular matrix1862 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.01867 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]  z1871 if j == i:1872 y = math.sqrt(z)1873 else:1874 y = z / fl[j][j]1875 fl[i][j] = y1876 #print_arr("ChoSol: fl", fl)1877 1878 # next, compute bl from fl and b1879 bl[0] = b[0] / fl[0][0]1880 for i in (1, 2):1881 z = 0.01882 for k in range(i):1883 z += fl[i][k] * bl[k]1884 #print "\t", i, k, z1885 bl[i] = (b[i]  z) / fl[i][i]1886 #print_vec("ChoSol: bl", bl)1887 1888 # last, compute beta from bl and fl1889 beta = np.empty((n))1890 beta[1] = bl[1] / fl[1][1]1891 for i in (1, 0):1892 z = 0.01893 for k in range(i+1, n):1894 z += fl[k][i] * beta[k]1895 #print "\t\t", i, k, 'z=', z1896 beta[i] = (bl[i]  z) / fl[i][i]1897 #print_vec("ChoSol: beta", beta)1898 1899 return beta1900 1901 def MaxEntMove(fSum, blank, chisq, chizer, c1, c2, s1, s2):1902 '''1903 move beta one step closer towards the solution1904 '''1905 a_lower, a_upper = 0., 1. # bracket "a"1906 cmin, beta = ChiNow (a_lower, c1, c2, s1, s2)1907 #print "MaxEntMove: cmin = %g" % cmin1908 if cmin*chisq > chizer:1909 ctarg = (1.0 + cmin)/21910 else:1911 ctarg = chizer/chisq1912 f_lower = cmin  ctarg1913 c_upper, beta = ChiNow (a_upper, c1, c2, s1, s2)1914 f_upper = c_upper  ctarg1915 1916 fx = 2*MOVE_PASSES # just to start off1917 loop = 11918 while abs(fx) >= MOVE_PASSES and loop <= MAX_MOVE_LOOPS:1919 a_new = (a_lower + a_upper) * 0.5 # search by bisection1920 c_new, beta = ChiNow (a_new, c1, c2, s1, s2)1921 fx = c_new  ctarg1922 # tighten the search range for the next pass1923 if f_lower*fx > 0:1924 a_lower, f_lower = a_new, fx1925 if f_upper*fx > 0:1926 a_upper, f_upper = a_new, fx1927 loop += 11928 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_DIRECTIONS1937 if (w > DISTANCE_LIMIT_FACTOR*fSum/blank): # invoke the distance penalty, SB eq. 171938 for k in range(m):1939 beta[k] *= math.sqrt (fSum/(blank*w))1940 chtarg = ctarg * chisq1941 return w, chtarg, loop, a_new, fx, beta1942 1943 #MaxEnt_SB starts here1944 1945 if image_to_data == None:1946 image_to_data = opus1947 if data_to_image == None:1948 data_to_image = tropus1949 n = len(base)1950 npt = len(datum)1951 1952 # Note that the order of subscripts for1953 # "xi" and "eta" has been reversed from1954 # the convention used in the FORTRAN version1955 # to enable parts of them to be passed as1956 # 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 base1967 1968 chizer, chtarg = npt*1.0, npt*1.01969 f = base * 1.0 # starting distribution is base1970 1971 fSum = sum(f) # find the sum of the fvector1972 z = (datum  image_to_data (f, G)) / sigma # standardized residuals, SB eq. 31973 chisq = sum(z*z) # Chi^2, SB eq. 41974 1975 for iter in range(IterMax):1976 ox = 2 * z / sigma # gradient of Chi^21977 1978 cgrad = data_to_image (ox, G) # cgrad[i] = del(C)/del(f[i]), SB eq. 81979 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. 221981 cnorm = math.sqrt(sum(f * cgrad*cgrad)) # ChiSqr term, SB eq. 221982 tnorm = sum(f * sgrad * cgrad) # norm for gradient term TEST1983 1984 a = 1.01985 b = 1.0 / cnorm1986 if iter == 0:1987 test = 0.0 # mismatch between entropy and ChiSquared gradients1988 else:1989 test = math.sqrt( ( 1.0  tnorm/(snorm*cnorm) )/2 ) # SB eq. 37?1990 a = 0.5 / (snorm * test)1991 b *= 0.5 / test1992 xi[0] = f * cgrad / cnorm1993 xi[1] = f * (a * sgrad  b * cgrad)1994 1995 eta[0] = image_to_data (xi[0], G); # image > data1996 eta[1] = image_to_data (xi[1], G); # image > data1997 ox = eta[1] / (sigma * sigma)1998 xi[2] = data_to_image (ox, G); # data > image1999 a = 1.0 / math.sqrt(sum(f * xi[2]*xi[2]))2000 xi[2] = f * xi[2] * a2001 eta[2] = image_to_data (xi[2], G) # image > data2002 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 technique2007 c1 = xi.dot(cgrad) / chisq # C_mu, SB eq. 242008 s1 = xi.dot(sgrad) # S_mu, SB eq. 242009 # 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) / chisq2015 s2[k][l] = sum(xi[k] * xi[l] / f) / blank2016 # print_arr("MaxEnt: c2", c2)2017 # print_arr("MaxEnt: s2", s2)2018 2019 # reflect across the body diagonal2020 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.02026 beta[2] = 0.02027 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 image2031 f += xi.transpose().dot(beta) # move the image towards the solution, SB eq. 252032 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_old2041 fSum = sum(f)2042 fChange = sum(df)2043 2044 # calculate the normalized entropy2045 S = sum((f/fSum) * np.log(f/fSum)) # normalized entropy, S&B eq. 12046 z = (datum  image_to_data (f, G)) / sigma # standardized residuals2047 chisq = sum(z*z) # report this ChiSq2048 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 = 02056 # 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 first2061 if (abs(chisq/chizer1.0) < CHI_SQR_LIMIT) and (test < TEST_LIMIT):2062 print ' Convergence achieved.'2063 return chisq,f,image_to_data(f, G) # solution FOUND returns here2064 print ' No convergence! Try increasing Error multiplier.'2065 return chisq,f,image_to_data(f, G) # no solution after IterMax iterations2066 2067 2068 ###############################################################################2069 #### IPG/TNNLS Routines2070 ###############################################################################2071 2072 def IPG(datum,sigma,G,Bins,Dbins,IterMax,Qvec=[],approach=0.8,Power=1,report=False):2073 ''' An implementation of the InteriorPoint Gradient method of2074 Michael Merritt & Yin Zhang, Technical Report TR0408, Dept. of Comp. and2075 Appl. Math., Rice Univ., Houston, Texas 77005, U.S.A. found on the web at2076 http://www.caam.rice.edu/caam/trs/2004/TR0408.pdf2077 Problem addressed: Total NonNegative Least Squares (TNNLS)2078 :param float datum[]:2079 :param float sigma[]:2080 :param float[][] G: transformation matrix2081 :param int IterMax:2082 :param float Qvec: data positions for Power = 042083 :param float approach: 0.8 default fitting parameter2084 :param int Power: 04 for Q^Power weighting, 1 to use input sigma2085 2086 '''2087 if Power < 0:2088 GmatE = G/sigma[:np.newaxis]2089 IntE = datum/sigma2090 pwr = 02091 QvecP = np.ones_like(datum)2092 else:2093 GmatE = G[:]2094 IntE = datum[:]2095 pwr = Power2096 QvecP = Qvec**pwr2097 Amat = GmatE*QvecP[:np.newaxis]2098 AAmat = np.inner(Amat,Amat)2099 Bvec = IntE*QvecP2100 Xw = np.ones_like(Bins)*1.e322101 calc = np.dot(G.T,Xw)2102 nIter = 02103 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*Qk2109 #Step 2 in M&Z paper:2110 alpSt = np.inner(Pk,Qk)/np.inner(Pk,np.inner(AAmat,Pk))2111 alpWv = Xw/Pk2112 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*Pk2115 Xw += shift2116 #done IPG; now results2117 nIter += 12118 calc = np.dot(G.T,Xw)2119 chisq = np.sum(((datumcalc)/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,calc2124 2125 ###############################################################################2126 #### SASD Utilities2127 ###############################################################################2128 2129 def SetScale(Data,refData):2130 Profile,Limits,Sample = Data2131 refProfile,refLimits,refSample = refData2132 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/sum2141 2142 def Bestimate(G,Rg,P):2143 return (G*P/Rg**P)*np.exp(scsp.gammaln(P/2))2144 2145 ###############################################################################2146 #### Size distribution2147 ###############################################################################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 radii2160 else:2161 Bins = np.linspace(data['Size']['MinDiam'],data['Size']['MaxDiam'],2162 data['Size']['Nbins']+1,True)/2. #make radii2163 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/Contrast2169 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 #### Modelling2195 ################################################################################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','SchulzZimm','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,values2244 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','SchulzZimm','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 M2272 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','SchulzZimm']: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**2rhoMat**22292 rBins,dBins,dist = MakeDiamDist(Type,parmDict[cid+'NumPoints'],parmDict[cid+'Cutoff'],distDict)2293 Gmat = G_matrix(Q,rBins,contrast,FFfxn,Volfxn,FFargs).T2294 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)))**32300 Guin = G*np.exp((Q*Rg)**2/3)2301 Porod = (B/Qstar**P)*np.exp((Q*Rgco)**2/3)2302 Ic += Guin+Porod2303 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 += Porod2307 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**2rhoMat**22316 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 Ic2323 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[:] = 02331 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/(IfinIbeglen(varyList)) #reduced chi^22339 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: ',IfinIbeg,' 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,covMatrix2350 except ValueError:2351 return False,0,0,0,0,02352 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[:] = 02373 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','SchulzZimm']: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**2rhoMat**22388 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).T2394 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)))**32403 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+Porod2406 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] += Porod2413 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**2rhoMat**22424 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. #pos2447 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([mode4.*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.,pos4.*distDict['StdDev']])2469 maxX = np.min([pos+4.*distDict['StdDev'],1.e5])2470 nP = 5002471 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 stuff2481 ################################################################################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.shape2495 print "%s[][] = (" % text2496 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 3column 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 columns2511 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, dI2515 print "MaxEnt_SB: "2516 test_data_file = os.path.join( 'testinp', 'test.sas')2517 rhosq = 100 # scattering contrast, 10^20 1/cm^42518 bkg = 0.1 # I = I  bkg2519 dMin, dMax, nRadii = 25, 9000, 402520 defaultDistLevel = 1.0e62521 IterMax = 402522 errFac = 1.052523 2524 r = np.logspace(math.log10(dMin), math.log10(dMax), nRadii)/22525 dr = r * (r[1]/r[0]  1) # step size2526 f_dr = np.ndarray((nRadii)) * 0 # volume fraction histogram2527 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 return2536 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.