Changeset 1249 for trunk/GSASIIsasd.py


Ignore:
Timestamp:
Mar 14, 2014 11:18:21 AM (8 years ago)
Author:
vondreele
Message:

include the 'sky' parameter in MaxEnt?; scale it by Scale & Contrast - makes it more universal.
Reorder the opus & tropus routines back to original definitions as per Pete J.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIsasd.py

    r1248 r1249  
    1313# $Id: GSASIIsasd.py 1186 2014-01-09 17:09:53Z vondreele $
    1414########### SVN repository information ###################
     15import os
    1516import sys
    1617import math
     
    283284    FF = FFfxn(q,r,args)
    284285    Vol = Volfxn(r,args)
    285     return 1.e-4*(contrast*Vol*FF**2)     #10^-20 vs 10^-24
     286    return 1.e-4*(contrast*Vol*FF**2).T     #10^-20 vs 10^-24
    286287   
    287288'''
     
    312313'''
    313314
    314 import os
    315 import sys
    316 import math
    317 import numpy
    318 
    319315class MaxEntException(Exception):
    320316    '''Any exception from this module'''
     
    345341    MOVE_PASSES       = 0.001                   # convergence test in routine: move
    346342
    347     def opus (data, G):
    348         '''
    349         opus: transform data-space -> solution-space:  [G] * data
     343    def tropus (data, G):
     344        '''
     345        tropus: transform data-space -> solution-space:  [G] * data
    350346       
    351347        default definition, caller can use this definition or provide an alternative
     
    357353        return G.dot(data)
    358354
    359     def tropus (image, G):
    360         '''
    361         tropus: transform solution-space -> data-space:  [G]^tr * image
     355    def opus (image, G):
     356        '''
     357        opus: transform solution-space -> data-space:  [G]^tr * image
    362358       
    363359        default definition, caller can use this definition or provide an alternative
     
    367363        :returns float[M]: calculated data, ndarray of shape (M)
    368364        '''
    369         return G.transpose().dot(image)
     365        return np.dot(G.T,image)    #G.transpose().dot(image)
    370366
    371367    def Dist(s2, beta):
     
    402398        '''
    403399        n = b.shape[0]
    404         fl = numpy.ndarray((n, n))*0
    405         bl = numpy.ndarray((n))*0
     400        fl = np.zeros((n,n))
     401        bl = np.zeros_like(b)
    406402       
    407403        #print_arr("ChoSol: a", a)
     
    443439   
    444440        # last, compute beta from bl and fl
    445         beta = numpy.ndarray((n))
     441        beta = np.empty((n))
    446442        beta[-1] = bl[-1] / fl[-1][-1]
    447443        for i in (1, 0):
     
    511507    # to enable parts of them to be passed as
    512508    # as vectors to "image_to_data" and "data_to_image".
    513     xi      = 0*numpy.ndarray((SEARCH_DIRECTIONS, n))
    514     eta     = 0*numpy.ndarray((SEARCH_DIRECTIONS, npt))
    515     beta    = 0*numpy.ndarray((SEARCH_DIRECTIONS))
    516     # s1      = 0*numpy.ndarray((SEARCH_DIRECTIONS))
    517     # c1      = 0*numpy.ndarray((SEARCH_DIRECTIONS))
    518     s2      = 0*numpy.ndarray((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
    519     c2      = 0*numpy.ndarray((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
     509    xi      = np.zeros((SEARCH_DIRECTIONS, n))
     510    eta     = np.zeros((SEARCH_DIRECTIONS, npt))
     511    beta    = np.zeros((SEARCH_DIRECTIONS))
     512    # s1      = np.zeros((SEARCH_DIRECTIONS))
     513    # c1      = np.zeros((SEARCH_DIRECTIONS))
     514    s2      = np.zeros((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
     515    c2      = np.zeros((SEARCH_DIRECTIONS, SEARCH_DIRECTIONS))
    520516
    521517    # TODO: replace blank (scalar) with base (vector)
     
    533529
    534530        cgrad = data_to_image (ox, G)              # cgrad[i] = del(C)/del(f[i]), SB eq. 8
    535         sgrad = -numpy.log(f/base) / (blank*math.exp (1.0))  # sgrad[i] = del(S)/del(f[i])
     531        sgrad = -np.log(f/base) / (blank*math.exp (1.0))  # sgrad[i] = del(S)/del(f[i])
    536532        snorm = math.sqrt(sum(f * sgrad*sgrad))    # entropy term, SB eq. 22
    537533        cnorm = math.sqrt(sum(f * cgrad*cgrad))    # ChiSqr term, SB eq. 22
     
    599595
    600596        # calculate the normalized entropy
    601         S = sum((f/fSum) * numpy.log(f/fSum))      # normalized entropy, S&B eq. 1
     597        S = sum((f/fSum) * np.log(f/fSum))      # normalized entropy, S&B eq. 1
    602598        z = (datum - image_to_data (f, G)) / sigma  # standardized residuals
    603599        chisq = sum(z*z)                            # report this ChiSq
     
    654650        M = len(buf)
    655651        buf = zip(*buf)         # transpose rows and columns
    656         q  = numpy.array(buf[0], dtype=numpy.float64)
    657         I  = numpy.array(buf[1], dtype=numpy.float64)
    658         dI = numpy.array(buf[2], dtype=numpy.float64)
     652        q  = np.array(buf[0], dtype=np.float64)
     653        I  = np.array(buf[1], dtype=np.float64)
     654        dI = np.array(buf[2], dtype=np.float64)
    659655        return q, I, dI
    660656    print "MaxEnt_SB: "
     
    667663    errFac = 1.05
    668664   
    669     r    = numpy.logspace(math.log10(dMin), math.log10(dMax), nRadii)/2
     665    r    = np.logspace(math.log10(dMin), math.log10(dMax), nRadii)/2
    670666    dr   = r * (r[1]/r[0] - 1)          # step size
    671     f_dr = numpy.ndarray((nRadii)) * 0  # volume fraction histogram
    672     b    = numpy.ndarray((nRadii)) * 0 + defaultDistLevel  # MaxEnt "sky background"
     667    f_dr = np.ndarray((nRadii)) * 0  # volume fraction histogram
     668    b    = np.ndarray((nRadii)) * 0 + defaultDistLevel  # MaxEnt "sky background"
    673669   
    674670    qVec, I, dI = readTextData(test_data_file)
     
    726722    Dbins = np.diff(Bins)
    727723    Bins = Bins[:-1]+Dbins/2.
    728     BinsBack = np.ones_like(Bins)*1.e-6
     724    Contrast = Sample['Contrast'][1]
     725    Scale = Sample['Scale'][0]
     726    Sky = 10**data['Size']['MaxEnt']['Sky']
     727    BinsBack = np.ones_like(Bins)*Sky*Scale/Contrast #How about *Scale/Contrast?
    729728    Back = data['Back']
    730729    Q,Io,wt,Ic,Ib = Profile[:5]
    731730    Qmin = Limits[1][0]
    732731    Qmax = Limits[1][1]
    733     Contrast = Sample['Contrast'][1]
    734732    wtFactor = ProfDict['wtFactor']
    735733    Ibeg = np.searchsorted(Q,Qmin)
     
    739737        Ic[Ibeg:Ifin] = Back[0]
    740738    Gmat = G_matrix(Q[Ibeg:Ifin],Bins,Contrast,shapes[Shape][0],shapes[Shape][1],args=Parms)
    741     chisq,BinMag,Ic[Ibeg:Ifin] = MaxEnt_SB(Sample['Scale'][0]*Io[Ibeg:Ifin]-Back[0],
    742         Sample['Scale'][0]/np.sqrt(wtFactor*wt[Ibeg:Ifin]),BinsBack,
     739    chisq,BinMag,Ic[Ibeg:Ifin] = MaxEnt_SB(Scale*Io[Ibeg:Ifin]-Back[0],
     740        Scale/np.sqrt(wtFactor*wt[Ibeg:Ifin]),BinsBack,
    743741        data['Size']['MaxEnt']['Niter'],Gmat,report=True)
    744742    print ' Final chi^2: %.3f'%(chisq)
Note: See TracChangeset for help on using the changeset viewer.