Changeset 940


Ignore:
Timestamp:
Jun 3, 2013 1:28:03 PM (9 years ago)
Author:
vondreele
Message:

conflicts fixed

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r939 r940  
    410410    return density,Volume/mass
    411411   
     412def getWave(Parms):
     413    try:
     414        return Parms['Lam'][1]
     415    except KeyError:
     416        return Parms['Lam1'][1]
     417   
     418################################################################################
     419##### distance, angle, planes, torsion stuff stuff
     420################################################################################
     421
    412422def getSyXYZ(XYZ,ops,SGData):
    413423    'Needs a doc string'
     
    9941004        out += ("e{:d}").format(valoff) # add an exponent, when needed
    9951005    return out
     1006
     1007################################################################################
     1008##### Fourier & charge flip stuff
     1009################################################################################
    9961010
    9971011def adjHKLmax(SGData,Hmax):
     
    14861500        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
    14871501    return XY
     1502   
     1503################################################################################
     1504##### MC/SA stuff
     1505################################################################################
     1506
     1507#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
     1508# Original Author: Travis Oliphant 2002
     1509# Bug-fixes in 2006 by Tim Leslie
     1510
     1511
     1512import numpy
     1513from numpy import asarray, tan, exp, ones, squeeze, sign, \
     1514     all, log, sqrt, pi, shape, array, minimum, where
     1515from numpy import random
     1516
     1517__all__ = ['anneal']
     1518
     1519_double_min = numpy.finfo(float).min
     1520_double_max = numpy.finfo(float).max
     1521class base_schedule(object):
     1522    def __init__(self):
     1523        self.dwell = 20
     1524        self.learn_rate = 0.5
     1525        self.lower = -10
     1526        self.upper = 10
     1527        self.Ninit = 50
     1528        self.accepted = 0
     1529        self.tests = 0
     1530        self.feval = 0
     1531        self.k = 0
     1532        self.T = None
     1533
     1534    def init(self, **options):
     1535        self.__dict__.update(options)
     1536        self.lower = asarray(self.lower)
     1537        self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
     1538        self.upper = asarray(self.upper)
     1539        self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
     1540        self.k = 0
     1541        self.accepted = 0
     1542        self.feval = 0
     1543        self.tests = 0
     1544
     1545    def getstart_temp(self, best_state):
     1546        """ Find a matching starting temperature and starting parameters vector
     1547        i.e. find x0 such that func(x0) = T0.
     1548
     1549        Parameters
     1550        ----------
     1551        best_state : _state
     1552            A _state object to store the function value and x0 found.
     1553
     1554        Returns
     1555        -------
     1556        x0 : array
     1557            The starting parameters vector.
     1558        """
     1559
     1560        assert(not self.dims is None)
     1561        lrange = self.lower
     1562        urange = self.upper
     1563        fmax = _double_min
     1564        fmin = _double_max
     1565        for _ in range(self.Ninit):
     1566            x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
     1567            fval = self.func(x0, *self.args)
     1568            self.feval += 1
     1569            if fval > fmax:
     1570                fmax = fval
     1571            if fval < fmin:
     1572                fmin = fval
     1573                best_state.cost = fval
     1574                best_state.x = array(x0)
     1575
     1576        self.T0 = (fmax-fmin)*1.5
     1577        return best_state.x
     1578
     1579    def accept_test(self, dE):
     1580        T = self.T
     1581        self.tests += 1
     1582        if dE < 0:
     1583            self.accepted += 1
     1584            return 1
     1585        p = exp(-dE*1.0/self.boltzmann/T)
     1586        if (p > random.uniform(0.0, 1.0)):
     1587            self.accepted += 1
     1588            return 1
     1589        return 0
     1590
     1591    def update_guess(self, x0):
     1592        pass
     1593
     1594    def update_temp(self, x0):
     1595        pass
     1596
     1597
     1598#  A schedule due to Lester Ingber
     1599class fast_sa(base_schedule):
     1600    def init(self, **options):
     1601        self.__dict__.update(options)
     1602        if self.m is None:
     1603            self.m = 1.0
     1604        if self.n is None:
     1605            self.n = 1.0
     1606        self.c = self.m * exp(-self.n * self.quench)
     1607
     1608    def update_guess(self, x0):
     1609        x0 = asarray(x0)
     1610        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
     1611        T = self.T
     1612        y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)
     1613        xc = y*(self.upper - self.lower)
     1614        xnew = x0 + xc
     1615        return xnew
     1616
     1617    def update_temp(self):
     1618        self.T = self.T0*exp(-self.c * self.k**(self.quench))
     1619        self.k += 1
     1620        return
     1621
     1622class cauchy_sa(base_schedule):
     1623    def update_guess(self, x0):
     1624        x0 = asarray(x0)
     1625        numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
     1626        xc = self.learn_rate * self.T * tan(numbers)
     1627        xnew = x0 + xc
     1628        return xnew
     1629
     1630    def update_temp(self):
     1631        self.T = self.T0/(1+self.k)
     1632        self.k += 1
     1633        return
     1634
     1635class boltzmann_sa(base_schedule):
     1636    def update_guess(self, x0):
     1637        std = minimum(sqrt(self.T)*ones(self.dims), (self.upper-self.lower)/3.0/self.learn_rate)
     1638        x0 = asarray(x0)
     1639        xc = squeeze(random.normal(0, 1.0, size=self.dims))
     1640
     1641        xnew = x0 + xc*std*self.learn_rate
     1642        return xnew
     1643
     1644    def update_temp(self):
     1645        self.k += 1
     1646        self.T = self.T0 / log(self.k+1.0)
     1647        return
     1648
     1649class log_sa(base_schedule):
     1650
     1651    def init(self,**options):
     1652        self.__dict__.update(options)
     1653       
     1654    def update_guess(self,x0):
     1655        x0 = np.asarray(x0)
     1656        u = np.squeeze(np.random.uniform(0.,1.,size=self.dims))
     1657        xnew = x0+u
     1658        return xnew
     1659       
     1660    def update_temp(self):
     1661        self.k += 1
     1662        self.T = self.T0*self.slope**k
     1663       
     1664class Tremayne_sa(base_schedule):   #needs fixing for two steps
     1665
     1666    def init(self,**options):
     1667        self.__dict__.update(options)
     1668
     1669    def update_guess(self,x0):
     1670        x0 = np.asarray(x0)
     1671        u = np.squeeze(np.random.uniform(0.,1.,size=self.dims))
     1672        xnew = x0+u
     1673        return xnew       
     1674   
     1675    def update_temp(self):
     1676        self.k += 1
     1677        self.T = self.T0*self.slope**k
     1678   
     1679class _state(object):
     1680    def __init__(self):
     1681        self.x = None
     1682        self.cost = None
     1683
     1684# TODO:
     1685#     allow for general annealing temperature profile
     1686#     in that case use update given by alpha and omega and
     1687#     variation of all previous updates and temperature?
     1688
     1689# Simulated annealing
     1690
     1691def anneal(func, x0, args=(), schedule='fast', full_output=0,
     1692           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
     1693           boltzmann=1.0, learn_rate=0.5, feps=1e-6, quench=1.0, m=1.0, n=1.0,
     1694           lower=-100, upper=100, dwell=50, slope=0.9):
     1695    """Minimize a function using simulated annealing.
     1696
     1697    Schedule is a schedule class implementing the annealing schedule.
     1698    Available ones are 'fast', 'cauchy', 'boltzmann'
     1699
     1700    Parameters
     1701    ----------
     1702    func : callable f(x, *args)
     1703        Function to be optimized.
     1704    x0 : ndarray
     1705        Initial guess.
     1706    args : tuple
     1707        Extra parameters to `func`.
     1708    schedule : base_schedule
     1709        Annealing schedule to use (a class).
     1710    full_output : bool
     1711        Whether to return optional outputs.
     1712    T0 : float
     1713        Initial Temperature (estimated as 1.2 times the largest
     1714        cost-function deviation over random points in the range).
     1715    Tf : float
     1716        Final goal temperature.
     1717    maxeval : int
     1718        Maximum function evaluations.
     1719    maxaccept : int
     1720        Maximum changes to accept.
     1721    maxiter : int
     1722        Maximum cooling iterations.
     1723    learn_rate : float
     1724        Scale constant for adjusting guesses.
     1725    boltzmann : float
     1726        Boltzmann constant in acceptance test
     1727        (increase for less stringent test at each temperature).
     1728    feps : float
     1729        Stopping relative error tolerance for the function value in
     1730        last four coolings.
     1731    quench, m, n : float
     1732        Parameters to alter fast_sa schedule.
     1733    lower, upper : float or ndarray
     1734        Lower and upper bounds on `x`.
     1735    dwell : int
     1736        The number of times to search the space at each temperature.
     1737    slope : float
     1738        Parameter for log schedule
     1739
     1740    Returns
     1741    -------
     1742    xmin : ndarray
     1743        Point giving smallest value found.
     1744    Jmin : float
     1745        Minimum value of function found.
     1746    T : float
     1747        Final temperature.
     1748    feval : int
     1749        Number of function evaluations.
     1750    iters : int
     1751        Number of cooling iterations.
     1752    accept : int
     1753        Number of tests accepted.
     1754    retval : int
     1755        Flag indicating stopping condition::
     1756
     1757                0 : Points no longer changing
     1758                1 : Cooled to final temperature
     1759                2 : Maximum function evaluations
     1760                3 : Maximum cooling iterations reached
     1761                4 : Maximum accepted query locations reached
     1762                5 : Final point not the minimum amongst encountered points
     1763
     1764    Notes
     1765    -----
     1766    Simulated annealing is a random algorithm which uses no derivative
     1767    information from the function being optimized. In practice it has
     1768    been more useful in discrete optimization than continuous
     1769    optimization, as there are usually better algorithms for continuous
     1770    optimization problems.
     1771
     1772    Some experimentation by trying the difference temperature
     1773    schedules and altering their parameters is likely required to
     1774    obtain good performance.
     1775
     1776    The randomness in the algorithm comes from random sampling in numpy.
     1777    To obtain the same results you can call numpy.random.seed with the
     1778    same seed immediately before calling scipy.optimize.anneal.
     1779
     1780    We give a brief description of how the three temperature schedules
     1781    generate new points and vary their temperature. Temperatures are
     1782    only updated with iterations in the outer loop. The inner loop is
     1783    over loop over xrange(dwell), and new points are generated for
     1784    every iteration in the inner loop. (Though whether the proposed
     1785    new points are accepted is probabilistic.)
     1786
     1787    For readability, let d denote the dimension of the inputs to func.
     1788    Also, let x_old denote the previous state, and k denote the
     1789    iteration number of the outer loop. All other variables not
     1790    defined below are input variables to scipy.optimize.anneal itself.
     1791
     1792    In the 'fast' schedule the updates are ::
     1793
     1794        u ~ Uniform(0, 1, size=d)
     1795        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
     1796        xc = y * (upper - lower)
     1797        x_new = x_old + xc
     1798
     1799        c = n * exp(-n * quench)
     1800        T_new = T0 * exp(-c * k**quench)
     1801
     1802
     1803    In the 'cauchy' schedule the updates are ::
     1804
     1805        u ~ Uniform(-pi/2, pi/2, size=d)
     1806        xc = learn_rate * T * tan(u)
     1807        x_new = x_old + xc
     1808
     1809        T_new = T0 / (1+k)
     1810
     1811    In the 'boltzmann' schedule the updates are ::
     1812
     1813        std = minimum( sqrt(T) * ones(d), (upper-lower) / (3*learn_rate) )
     1814        y ~ Normal(0, std, size=d)
     1815        x_new = x_old + learn_rate * y
     1816
     1817        T_new = T0 / log(1+k)
     1818
     1819    """
     1820    x0 = asarray(x0)
     1821    lower = asarray(lower)
     1822    upper = asarray(upper)
     1823
     1824    schedule = eval(schedule+'_sa()')
     1825    #   initialize the schedule
     1826    schedule.init(dims=shape(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
     1827                  learn_rate=learn_rate, lower=lower, upper=upper,
     1828                  m=m, n=n, quench=quench, dwell=dwell)
     1829
     1830    current_state, last_state, best_state = _state(), _state(), _state()
     1831    if T0 is None:
     1832        x0 = schedule.getstart_temp(best_state)
     1833    else:
     1834        best_state.x = None
     1835        best_state.cost = numpy.Inf
     1836
     1837    last_state.x = asarray(x0).copy()
     1838    fval = func(x0,*args)
     1839    schedule.feval += 1
     1840    last_state.cost = fval
     1841    if last_state.cost < best_state.cost:
     1842        best_state.cost = fval
     1843        best_state.x = asarray(x0).copy()
     1844    schedule.T = schedule.T0
     1845    fqueue = [100, 300, 500, 700]
     1846    iters = 0
     1847    while 1:
     1848        for n in xrange(dwell):
     1849            current_state.x = schedule.update_guess(last_state.x)
     1850            current_state.cost = func(current_state.x,*args)
     1851            schedule.feval += 1
     1852
     1853            dE = current_state.cost - last_state.cost
     1854            if schedule.accept_test(dE):
     1855                last_state.x = current_state.x.copy()
     1856                last_state.cost = current_state.cost
     1857                if last_state.cost < best_state.cost:
     1858                    best_state.x = last_state.x.copy()
     1859                    best_state.cost = last_state.cost
     1860        schedule.update_temp()
     1861        iters += 1
     1862        # Stopping conditions
     1863        # 0) last saved values of f from each cooling step
     1864        #     are all very similar (effectively cooled)
     1865        # 1) Tf is set and we are below it
     1866        # 2) maxeval is set and we are past it
     1867        # 3) maxiter is set and we are past it
     1868        # 4) maxaccept is set and we are past it
     1869
     1870        fqueue.append(squeeze(last_state.cost))
     1871        fqueue.pop(0)
     1872        af = asarray(fqueue)*1.0
     1873        if all(abs((af-af[0])/af[0]) < feps):
     1874            retval = 0
     1875            if abs(af[-1]-best_state.cost) > feps*10:
     1876                retval = 5
     1877                print "Warning: Cooled to %f at %s but this is not" \
     1878                      % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
     1879                      + " the smallest point found."
     1880            break
     1881        if (Tf is not None) and (schedule.T < Tf):
     1882            retval = 1
     1883            break
     1884        if (maxeval is not None) and (schedule.feval > maxeval):
     1885            retval = 2
     1886            break
     1887        if (iters > maxiter):
     1888            print "Warning: Maximum number of iterations exceeded."
     1889            retval = 3
     1890            break
     1891        if (maxaccept is not None) and (schedule.accepted > maxaccept):
     1892            retval = 4
     1893            break
     1894
     1895    if full_output:
     1896        return best_state.x, best_state.cost, schedule.T, \
     1897               schedule.feval, iters, schedule.accepted, retval
     1898    else:
     1899        return best_state.x, retval
     1900
    14881901
    14891902def mcsaSearch(data,reflType,reflData,covData,pgbar):
    1490     'Needs a doc string'       
     1903    gamFW = lambda s,g: math.exp(math.log(s**5+2.69269*s**4*g+2.42843*s**3*g**2+4.47163*s**2*g**3+0.07842*s*g**4+g**5)/5.)
     1904   
     1905    def getMDparms(item,pfx,parmDict,varyList):
     1906        parmDict[pfx+'MDaxis'] = item['axis']
     1907        parmDict[pfx+'MDval'] = item['Coef'][0]
     1908        if item['Coef'][1]:
     1909            varyList += [pfx+'MDval',]
     1910            limits = item['Coef'][2]
     1911            lower.append(limits[0])
     1912            upper.append(limits[1])
     1913           
     1914    def getAtomparms(item,pfx,parmDict,varyList):
     1915        parmDict[pfx+'Atype'] = item['atType']
     1916        pstr = ['x','y','z']
     1917        for i in range(3):
     1918            name = pfx+pstr[i]
     1919            parmDict[name] = item['Pos'][0][i]
     1920            if item['Pos'][1][i]:
     1921                varyList += [name,]
     1922                limits = item['Pos'][2][i]
     1923                lower.append(limits[0])
     1924                upper.append(limits[1])
     1925           
     1926    def getRBparms(item,pfx,parmDict,varyList):
     1927        parmDict['MolCent'] = item['MolCent']
     1928        parmDict['RBId'] = item['RBId']
     1929        pstr = ['Px','Py','Pz']
     1930        ostr = ['Qa','Qi','Qj','Qk']
     1931        for i in range(3):
     1932            name = pfx+pstr[i]
     1933            parmDict[name] = item['Pos'][0][i]
     1934            if item['Pos'][1][i]:
     1935                varyList += [name,]
     1936                limits = item['Pos'][2][i]
     1937                lower.append(limits[0])
     1938                upper.append(limits[1])
     1939        for i in range(4):
     1940            name = pfx+ostr[i]
     1941            parmDict[name] = item['Ori'][0][i]
     1942            if item['Ovar'] == 'AV' and i:
     1943                varyList += [name,]
     1944                limits = item['Ori'][2][i]
     1945                lower.append(limits[0])
     1946                upper.append(limits[1])
     1947            elif item['Ovar'] == 'A' and not i:
     1948                varyList += [name,]
     1949                limits = item['Ori'][2][i]
     1950                lower.append(limits[0])
     1951                upper.append(limits[1])
     1952        if 'Tor' in item:      #'Tor' not there for 'Vector' RBs
     1953            for i in range(len(item['Tor'][0])):
     1954                name = pfx+'Tor'+str(i)
     1955                parmDict[name] = item['Tor'][0][i]
     1956                if item['Tor'][1][i]:
     1957                    varyList += [name,]
     1958                    limits = item['Tor'][2][i]
     1959                    lower.append(limits[0])
     1960                    upper.append(limits[1])
     1961
     1962    def getFFvalues(FFtables,SQ):
     1963        FFvals = {}
     1964        for El in FFtables:
     1965            FFvals[El] = G2el.ScatFac(FFtables[El],SQ)[0]
     1966        return FFvals
     1967       
     1968    def getBLvalues(BLtables):
     1969        BLvals = {}
     1970        for El in BLtables:
     1971            BLvals[El] = BLtables[El][1][1]
     1972        return BLvals
     1973           
     1974#    def mcsaSfCalc(refList,G,SGData,calcControls,parmDict):
     1975#        ''' Compute structure factors for all h,k,l for phase
     1976#        input:
     1977#            refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
     1978#            G:      reciprocal metric tensor
     1979#            SGData: space group info. dictionary output from SpcGroup
     1980#            calcControls:
     1981#            ParmDict:
     1982#        puts result F^2 in each ref[8] in refList
     1983#        '''       
     1984#        twopi = 2.0*np.pi
     1985#        twopisq = 2.0*np.pi**2
     1986#        ast = np.sqrt(np.diag(G))
     1987#        Mast = twopisq*np.multiply.outer(ast,ast)
     1988#        FFtables = calcControls['FFtables']
     1989#        BLtables = calcControls['BLtables']
     1990#        Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     1991#        FF = np.zeros(len(Tdata))
     1992#        if 'N' in parmDict[hfx+'Type']:
     1993#            FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
     1994#        else:
     1995#            FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
     1996#            FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
     1997#        maxPos = len(SGData['SGOps'])
     1998#        Uij = np.array(G2lat.U6toUij(Uijdata))
     1999#        bij = Mast*Uij.T
     2000#        for refl in refList:
     2001#            fbs = np.array([0,0])
     2002#            H = refl[:3]
     2003#            SQ = 1./(2.*refl[4])**2
     2004#            SQfactor = 4.0*SQ*twopisq
     2005#            Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
     2006#            if not len(refl[-1]):                #no form factors
     2007#                if 'N' in parmDict[hfx+'Type']:
     2008#                    refl[-1] = getBLvalues(BLtables)
     2009#                else:       #'X'
     2010#                    refl[-1] = getFFvalues(FFtables,SQ)
     2011#            for i,El in enumerate(Tdata):
     2012#                FF[i] = refl[-1][El]           
     2013#            Uniq = refl[11]
     2014#            phi = refl[12]
     2015#            phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
     2016#            sinp = np.sin(phase)
     2017#            cosp = np.cos(phase)
     2018#            occ = Mdata*Fdata/len(Uniq)
     2019#            biso = -SQfactor*Uisodata
     2020#            Tiso = np.where(biso<1.,np.exp(biso),1.0)
     2021#            HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
     2022#            Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
     2023#            Tcorr = Tiso*Tuij
     2024#            fa = np.array([(FF+FP-Bab)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr])
     2025#            fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
     2026#            if not SGData['SGInv']:
     2027#                fb = np.array([(FF+FP-Bab)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
     2028#                fbs = np.sum(np.sum(fb,axis=1),axis=1)
     2029#            fasq = fas**2
     2030#            fbsq = fbs**2        #imaginary
     2031#            refl[9] = np.sum(fasq)+np.sum(fbsq)
     2032#            refl[10] = atan2d(fbs[0],fas[0])
     2033   
     2034    sq8ln2 = np.sqrt(8*np.log(2))
     2035    sq2pi = np.sqrt(2*np.pi)
     2036    sq4pi = np.sqrt(4*np.pi)
    14912037    generalData = data['General']
     2038    pId = data['pId']
     2039    pfx = str(pId)+'::'
     2040    Atoms = data['Atoms']                       #if any
    14922041    mcsaControls = generalData['MCSA controls']
    1493     reflName = mcsaData['Data source']
     2042    reflName = mcsaControls['Data source']
    14942043    phaseName = generalData['Name']
    1495     MCSAObjs = data['MCSAModels']
     2044    MCSAObjs = data['MCSA']['Models']               #list of MCSA models
     2045    parmDict = {}
     2046    varyList = []
     2047    upper = []
     2048    lower = []
     2049    for i,item in enumerate(MCSAObjs):
     2050        mfx = str(i)+':'
     2051        if item['Type'] == 'MD':
     2052            getMDparms(item,mfx,parmDict,varyList)
     2053        elif item['Type'] == 'Atom':
     2054            getAtomparms(item,mfx,parmDict,varyList)
     2055        elif item['Type'] in ['Residue','Vector']:
     2056            getRBparms(item,mfx,parmDict,varyList)
     2057    if 'PWDR' in reflName:
     2058        refs = []
     2059        for ref in reflData:
     2060            h,k,l,m,d,pos,sig,gam,Fsq = ref[:9]
     2061            if d >= mcsaControls['dmin']:
     2062                FWHM = gamFW(sig,gam)/2.35482        #sqrt(8ln2) --> sig from FWHM
     2063                Fsq *= m
     2064    elif 'Pawley' in reflName:
     2065        refs = []
     2066        covMatrix = covData['covMatrix']
     2067        varyList = covData['varyList']
     2068        for iref,refI in enumerate(reflData):
     2069            h,k,l,m,d,v,f,s = refI
     2070            if d >= mcsaControls['dmin'] and v:       #skip unrefined ones
     2071                nameI = pfx+'PWLref:'+str(iref)
     2072                if nameI in covData['varyList']:
     2073                    refs.append([h,k,l,f*m,0.,0.])
     2074        nRef = len(refs)
     2075        print nRef       
     2076        covTerms = np.zeros((nRef,nRef))
     2077        print covTerms.shape
     2078       
     2079        Iref = 0
     2080        for iref,refI in enumerate(reflData):
     2081            if refI[4] >= mcsaControls['dmin'] and refI[5]:       #skip unrefined ones
     2082                nameI = pfx+'PWLref:'+str(iref)
     2083                if nameI in covData['varyList']:
     2084                    Iindx = varyList.index(nameI)
     2085                    covTerms[Iref][Iref] = covMatrix[Iindx][Iindx]
     2086                    Jref = 0
     2087                    for jref,refJ in enumerate(reflData[:iref]):
     2088                        if refJ[5]:
     2089                            nameJ = pfx+'PWLref:'+str(jref)
     2090                            try:
     2091                                covTerms[Iref][Jref] = covMatrix[varyList.index(nameI)][varyList.index(nameJ)]
     2092                            except ValueError:
     2093                                covTerms[Iref][Jref] = covTerms[Iref][Jref-1]
     2094                            Jref += 1
     2095                Iref += 1
     2096        print covTerms
     2097                       
     2098                   
     2099#    covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
     2100#        'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
     2101    elif 'HKLF' in reflName:
     2102        refs = []
     2103        for ref in reflData:
     2104            h,k,l,m,d,Fsq = ref[:5],ref[6]
     2105            if d >= mcsaControls['dmin']:
     2106                Fsq *= m
     2107                refs.append([h,k,l,m,Fsq])
     2108        rcov = np.identity(len(refs))   
     2109                       
     2110       
    14962111           
    1497 #             = {'':'','Annealing':[50.0,0.001,0.90,1000],
    1498 #            'dmin':2.0,'Algolrithm':'Normal','Jump coeff':[0.95,0.5]} #'Normal','Random jump','Tremayne jump'
     2112#    generalData['MCSA controls'] = {'Data source':'','Annealing':[50.,0.001,50,1.e-6],
     2113#    'dmin':2.0,'Algorithm':'fast','Jump coeff':[0.95,0.5],'nRuns':1,'boltzmann':1.0,
     2114#    'fast parms':[1.0,1.0,1.0],'log slope':0.9}
    14992115
    15002116    return {}
    15012117
    15022118       
    1503 def getWave(Parms):
    1504     'Needs a doc string'       
    1505     try:
    1506         return Parms['Lam'][1]
    1507     except KeyError:
    1508         return Parms['Lam1'][1]
    1509    
     2119################################################################################
     2120##### Quaternion stuff
     2121################################################################################
     2122
    15102123def prodQQ(QA,QB):
    15112124    ''' Grassman quaternion product
     
    16632276            D *= -1.
    16642277    return Q,D
     2278   
     2279if __name__ == "__main__":
     2280    from numpy import cos
     2281    # minimum expected at ~-0.195
     2282    func = lambda x: cos(14.5*x-0.3) + (x+0.2)*x
     2283    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='cauchy')
     2284    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='fast')
     2285    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='boltzmann')
     2286
     2287    # minimum expected at ~[-0.195, -0.1]
     2288    func = lambda x: cos(14.5*x[0]-0.3) + (x[1]+0.2)*x[1] + (x[0]+0.2)*x[0]
     2289    print anneal(func,[1.0, 1.0],full_output=1,upper=[3.0, 3.0],lower=[-3.0, -3.0],feps=1e-4,maxiter=2000,schedule='cauchy')
     2290    print anneal(func,[1.0, 1.0],full_output=1,upper=[3.0, 3.0],lower=[-3.0, -3.0],feps=1e-4,maxiter=2000,schedule='fast')
     2291    print anneal(func,[1.0, 1.0],full_output=1,upper=[3.0, 3.0],lower=[-3.0, -3.0],feps=1e-4,maxiter=2000,schedule='boltzmann')
  • trunk/GSASIIphsGUI.py

    r936 r940  
    120120        if 'Pawley neg wt' not in generalData:
    121121            generalData['Pawley neg wt'] = 0.0
    122         if 'MCSA controls' not in generalData:
    123             generalData['MCSA controls'] = {'Data source':'','Annealing':[50.0,0.001,0.90,1000],
    124             'dmin':2.0,'Algolrithm':'Random jump','Jump coeff':[0.95,0.5],'nRuns':1} #'Random jump','Tremayne jump'
     122        if 'Algolrithm' in generalData.get('MCSA controls',{}) or \
     123            'MCSA controls' not in generalData:
     124            generalData['MCSA controls'] = {'Data source':'','Annealing':[50.,0.001,50,1.e-6],
     125            'dmin':2.0,'Algorithm':'fast','Jump coeff':[0.95,0.5],'nRuns':1,'boltzmann':1.0,
     126            'fast parms':[1.0,1.0,1.0],'log slope':0.9}
    125127# end of patches
    126128        generalData['NoAtoms'] = {}
     
    656658           
    657659            def OnAlist(event):
    658                 MCSA['Algolrithm'] = Alist.GetValue()
    659                 if 'Tremayne' in MCSA['Algolrithm']:
    660                     wx.CallAfter(UpdateGeneral)
     660                MCSA['Algorithm'] = Alist.GetValue()
     661                wx.CallAfter(UpdateGeneral)
     662               
     663            def OnSlope(event):
     664                try:
     665                    val = float(slope.GetValue())
     666                    if .25 <= val < 1.0:
     667                        MCSA['log slope'] = val
     668                except ValueError:
     669                    pass
     670                Obj.SetValue("%.3f"%(MCSA['log slope']))          #reset in case of error               
    661671           
    662672            def OnAjump(event):
    663673                Obj = event.GetEventObject()
    664                 ind = Indx[Obj.GetId()]
     674                name,ind = Indx[Obj.GetId()]
    665675                try:
    666676                    val = float(Obj.GetValue())
    667677                    if .0 <= val <= 1.0:
    668                         MCSA['Jump coeff'][ind] = val
     678                        MCSA[name][ind] = val
    669679                except ValueError:
    670680                    pass
    671                 Obj.SetValue("%.3f"%(MCSA['Jump coeff'][ind]))          #reset in case of error
     681                Obj.SetValue("%.3f"%(MCSA[name][ind]))
    672682           
    673683            def OnAnneal(event):
    674684                Obj = event.GetEventObject()
    675685                ind,fmt = Indx[Obj.GetId()]
    676                 if ind < 3:
     686                if ind == 2:        #No. trials
     687                    try:
     688                        val = int(Obj.GetValue())
     689                        if 1 <= val:
     690                            MCSA['Annealing'][ind] = val
     691                    except ValueError:
     692                        pass
     693                else:
    677694                    try:
    678695                        val = float(Obj.GetValue())
     
    680697                            MCSA['Annealing'][ind] = val
    681698                    except ValueError:
    682                         pass
    683                 else:
    684                     try:
    685                         val = int(Obj.GetValue())
    686                         if 1 <= val:
    687                             MCSA['Annealing'][3] = val
    688                     except ValueError:
    689                         pass
    690                    
    691                 Obj.SetValue(fmt%(MCSA['Annealing'][ind]))          #reset in case of error
     699                        pass                   
     700                Obj.SetValue(fmt%(MCSA['Annealing'][ind]))
    692701                       
     702#            MCSA = {'boltzmann':1.0,
    693703            refList = []
    694704            if len(data['Pawley ref']):
     
    712722            mcsaSizer.Add((5,5),)
    713723            line2Sizer = wx.BoxSizer(wx.HORIZONTAL)
    714             Rchoice = ['1','2','3','5','10','15','20']
     724            Rchoice = ['1','2','3','5','10','15','20','50','100','200','500']
    715725            line2Sizer.Add(wx.StaticText(General,label=' No. MC/SA runs: '),0,wx.ALIGN_CENTER_VERTICAL)
    716726            noRuns = wx.ComboBox(General,-1,value=str(MCSA.get('nRuns',1)),choices=Rchoice,
     
    718728            noRuns.Bind(wx.EVT_COMBOBOX,OnNoRuns)
    719729            line2Sizer.Add(noRuns,0,wx.ALIGN_CENTER_VERTICAL)
    720             Achoice = ['Random jump','Tremayne jump']
    721             line2Sizer.Add(wx.StaticText(General,label=' MC/SA algorithm: '),0,wx.ALIGN_CENTER_VERTICAL)
    722             Alist = wx.ComboBox(General,-1,value=MCSA['Algolrithm'],choices=Achoice,
     730            Achoice = ['log','fast','cauchy','boltzmann','Tremayne']
     731            line2Sizer.Add(wx.StaticText(General,label=' MC/SA schedule: '),0,wx.ALIGN_CENTER_VERTICAL)
     732            Alist = wx.ComboBox(General,-1,value=MCSA['Algorithm'],choices=Achoice,
    723733                style=wx.CB_READONLY|wx.CB_DROPDOWN)
    724734            Alist.Bind(wx.EVT_COMBOBOX,OnAlist)
    725735            line2Sizer.Add(Alist,0,wx.ALIGN_CENTER_VERTICAL)
    726             if 'Tremayne' in MCSA['Algolrithm']:
    727                 for i,name in enumerate([' A-jump: ',' B-jump: ']):
     736            if MCSA['Algorithm'] in ['Tremayne','fast','boltzmann','cauchy']:
     737                Names = [' A-jump: ',' B-jump: ']
     738                parms = 'Jump coeff'
     739                if MCSA['Algorithm'] in ['boltzmann','cauchy']:
     740                    Names = [' A-jump: ']
     741                elif 'fast' in MCSA['Algorithm']:
     742                    Names = [' quench: ',' m-factor: ',' n-factor: ']
     743                    parms = 'fast parms'
     744                for i,name in enumerate(Names):
    728745                    line2Sizer.Add(wx.StaticText(General,label=name),0,wx.ALIGN_CENTER_VERTICAL)
    729                     Ajump =  wx.TextCtrl(General,-1,value='%.3f'%(MCSA['Jump coeff'][i]),style=wx.TE_PROCESS_ENTER)
     746                    Ajump =  wx.TextCtrl(General,-1,value='%.3f'%(MCSA[parms][i]),style=wx.TE_PROCESS_ENTER)
    730747                    Ajump.Bind(wx.EVT_TEXT_ENTER,OnAjump)       
    731748                    Ajump.Bind(wx.EVT_KILL_FOCUS,OnAjump)
    732                     Indx[Ajump.GetId()] = i
     749                    Indx[Ajump.GetId()] = [parms,i]
    733750                    line2Sizer.Add(Ajump,0,wx.ALIGN_CENTER_VERTICAL)
     751            elif 'log' in MCSA['Algorithm']:
     752                line2Sizer.Add(wx.StaticText(General,label=' slope: '),0,wx.ALIGN_CENTER_VERTICAL)
     753                slope =  wx.TextCtrl(General,-1,value='%.3f'%(MCSA['log slope']),style=wx.TE_PROCESS_ENTER)
     754                slope.Bind(wx.EVT_TEXT_ENTER,OnSlope)       
     755                slope.Bind(wx.EVT_KILL_FOCUS,OnSlope)
     756                line2Sizer.Add(slope,0,wx.ALIGN_CENTER_VERTICAL)
    734757            mcsaSizer.Add(line2Sizer)
    735758            mcsaSizer.Add((5,5),)
    736759            line3Sizer = wx.BoxSizer(wx.HORIZONTAL)
    737760            line3Sizer.Add(wx.StaticText(General,label=' Annealing schedule: '),0,wx.ALIGN_CENTER_VERTICAL)
    738             names = [' Start temp: ',' Final temp: ',' Slope: ',' No. trials: ']
    739             fmts = ['%.1f','%.5f','%.2f','%d']
     761            names = [' Start temp: ',' Final temp: ',' No. trials: ',' feps: ']
     762            fmts = ['%.1f','%.5f','%d','%.2g']
    740763            for i,[name,fmt] in enumerate(zip(names,fmts)):
    741764                line3Sizer.Add(wx.StaticText(General,label=name),0,wx.ALIGN_CENTER_VERTICAL)
Note: See TracChangeset for help on using the changeset viewer.