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