Changeset 940
- Timestamp:
- Jun 3, 2013 1:28:03 PM (10 years ago)
- Location:
- trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIImath.py
r939 r940 410 410 return density,Volume/mass 411 411 412 def 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 412 422 def getSyXYZ(XYZ,ops,SGData): 413 423 'Needs a doc string' … … 994 1004 out += ("e{:d}").format(valoff) # add an exponent, when needed 995 1005 return out 1006 1007 ################################################################################ 1008 ##### Fourier & charge flip stuff 1009 ################################################################################ 996 1010 997 1011 def adjHKLmax(SGData,Hmax): … … 1486 1500 XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0] 1487 1501 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 1512 import numpy 1513 from numpy import asarray, tan, exp, ones, squeeze, sign, \ 1514 all, log, sqrt, pi, shape, array, minimum, where 1515 from numpy import random 1516 1517 __all__ = ['anneal'] 1518 1519 _double_min = numpy.finfo(float).min 1520 _double_max = numpy.finfo(float).max 1521 class 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 1599 class 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 1622 class 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 1635 class 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 1649 class 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 1664 class 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 1679 class _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 1691 def 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 1488 1901 1489 1902 def 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) 1491 2037 generalData = data['General'] 2038 pId = data['pId'] 2039 pfx = str(pId)+'::' 2040 Atoms = data['Atoms'] #if any 1492 2041 mcsaControls = generalData['MCSA controls'] 1493 reflName = mcsa Data['Data source']2042 reflName = mcsaControls['Data source'] 1494 2043 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 1496 2111 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} 1499 2115 1500 2116 return {} 1501 2117 1502 2118 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 1510 2123 def prodQQ(QA,QB): 1511 2124 ''' Grassman quaternion product … … 1663 2276 D *= -1. 1664 2277 return Q,D 2278 2279 if __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 120 120 if 'Pawley neg wt' not in generalData: 121 121 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} 125 127 # end of patches 126 128 generalData['NoAtoms'] = {} … … 656 658 657 659 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 661 671 662 672 def OnAjump(event): 663 673 Obj = event.GetEventObject() 664 ind = Indx[Obj.GetId()]674 name,ind = Indx[Obj.GetId()] 665 675 try: 666 676 val = float(Obj.GetValue()) 667 677 if .0 <= val <= 1.0: 668 MCSA[ 'Jump coeff'][ind] = val678 MCSA[name][ind] = val 669 679 except ValueError: 670 680 pass 671 Obj.SetValue("%.3f"%(MCSA[ 'Jump coeff'][ind])) #reset in case of error681 Obj.SetValue("%.3f"%(MCSA[name][ind])) 672 682 673 683 def OnAnneal(event): 674 684 Obj = event.GetEventObject() 675 685 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: 677 694 try: 678 695 val = float(Obj.GetValue()) … … 680 697 MCSA['Annealing'][ind] = val 681 698 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])) 692 701 702 # MCSA = {'boltzmann':1.0, 693 703 refList = [] 694 704 if len(data['Pawley ref']): … … 712 722 mcsaSizer.Add((5,5),) 713 723 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'] 715 725 line2Sizer.Add(wx.StaticText(General,label=' No. MC/SA runs: '),0,wx.ALIGN_CENTER_VERTICAL) 716 726 noRuns = wx.ComboBox(General,-1,value=str(MCSA.get('nRuns',1)),choices=Rchoice, … … 718 728 noRuns.Bind(wx.EVT_COMBOBOX,OnNoRuns) 719 729 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['Algo lrithm'],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, 723 733 style=wx.CB_READONLY|wx.CB_DROPDOWN) 724 734 Alist.Bind(wx.EVT_COMBOBOX,OnAlist) 725 735 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): 728 745 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) 730 747 Ajump.Bind(wx.EVT_TEXT_ENTER,OnAjump) 731 748 Ajump.Bind(wx.EVT_KILL_FOCUS,OnAjump) 732 Indx[Ajump.GetId()] = i749 Indx[Ajump.GetId()] = [parms,i] 733 750 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) 734 757 mcsaSizer.Add(line2Sizer) 735 758 mcsaSizer.Add((5,5),) 736 759 line3Sizer = wx.BoxSizer(wx.HORIZONTAL) 737 760 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'] 740 763 for i,[name,fmt] in enumerate(zip(names,fmts)): 741 764 line3Sizer.Add(wx.StaticText(General,label=name),0,wx.ALIGN_CENTER_VERTICAL)
Note: See TracChangeset
for help on using the changeset viewer.