Changeset 1301 for trunk/GSASIIsasd.py


Ignore:
Timestamp:
Apr 25, 2014 12:51:35 PM (8 years ago)
Author:
vondreele
Message:
 
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 ###################
     16import os
     17import sys
     18import math
     19import time
     20
     21import numpy as np
     22import scipy as sp
     23import numpy.linalg as nl
     24from numpy.fft import ifft, fft, fftshift
     25import scipy.special as scsp
     26import scipy.interpolate as si
     27import scipy.stats as st
     28import scipy.optimize as so
     29#import pdb
     30
     31import GSASIIpath
     32GSASIIpath.SetVersionNumber("$Revision$")
     33import GSASIIlattice as G2lat
     34import GSASIIspc as G2spc
     35import GSASIIElem as G2elem
     36import GSASIIgrid as G2gd
     37import GSASIIIO as G2IO
     38import GSASIImath as G2mth
     39import GSASIIpwd as G2pwd
     40
     41# trig functions in degrees
     42sind = lambda x: math.sin(x*math.pi/180.)
     43asind = lambda x: 180.*math.asin(x)/math.pi
     44tand = lambda x: math.tan(x*math.pi/180.)
     45atand = lambda x: 180.*math.atan(x)/math.pi
     46atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
     47cosd = lambda x: math.cos(x*math.pi/180.)
     48acosd = lambda x: 180.*math.acos(x)/math.pi
     49rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
     50#numpy versions
     51npsind = lambda x: np.sin(x*np.pi/180.)
     52npasind = lambda x: 180.*np.arcsin(x)/math.pi
     53npcosd = lambda x: np.cos(x*math.pi/180.)
     54npacosd = lambda x: 180.*np.arccos(x)/math.pi
     55nptand = lambda x: np.tan(x*math.pi/180.)
     56npatand = lambda x: 180.*np.arctan(x)/np.pi
     57npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
     58npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave
     59npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)
     60   
     61###############################################################################
     62#### Particle form factors
     63###############################################################################
     64
     65def 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   
     75def 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           
     92def 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   
     111def 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   
     121def 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   
     131def 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   
     143def 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   
     162def 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   
     172def 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   
     192def 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
     221def 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
     230def 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   
     240def 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   
     250def 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   
     260def 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   
     270def 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   
     279def 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   
     289def 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   
     299def 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   
     309def 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
     323def 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   
     335def 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   
     346def 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   
     360def 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           
     377def 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   
     390def 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   
     402def 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   
     418def 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
     443def 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'''
     458sbmaxent
     459
     460Entropy maximization routine as described in the article
     461J 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
     468References:
     469
     4701. J Skilling and RK Bryan; MON NOT R ASTR SOC 211 (1984) 111 - 124.
     4712. 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.
     4763. ID Culverwell and GP Clarke; Ibid. 87 - 96.
     4774. JA Potton, GK Daniell, & BD Rainford,
     478   J APPL CRYST 21 (1988) 663 - 668.
     4795. JA Potton, GJ Daniell, & BD Rainford,
     480   J APPL CRYST 21 (1988) 891 - 897.
     481
     482'''
     483
     484class MaxEntException(Exception):
     485    '''Any exception from this module'''
     486    pass
     487
     488def 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
     793def 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
     850def 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   
     863def Bestimate(G,Rg,P):
     864    return (G*P/Rg**P)*np.exp(scsp.gammaln(P/2))
     865   
     866###############################################################################
     867#### Size distribution
     868###############################################################################
     869
     870def 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
     918def 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   
     1086def 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   
     1170def 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
     1216def 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
     1225def 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
     1236def 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
     1274def tests():
     1275    test_MaxEnt_SB(report=True)
     1276
     1277if __name__ == '__main__':
     1278    tests()
     1279
     1280=======
    11281#/usr/bin/env python
    21282# -*- coding: utf-8 -*-
     
    12652545    tests()
    12662546
     2547>>>>>>> .r1299
Note: See TracChangeset for help on using the changeset viewer.