Changeset 4025
- Timestamp:
- Jun 13, 2019 3:56:48 PM (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIImath.py
r4016 r4025 29 29 import GSASIIspc as G2spc 30 30 import GSASIIpwd as G2pwd 31 import GSASIIfiles as G2fil 31 32 import numpy.fft as fft 32 33 import scipy.optimize as so … … 157 158 nfev = 0 158 159 if Print: 159 print(' Hessian Levenberg-Marquardt SVD refinement on %d variables:'%(n))160 G2fil.G2Print(' Hessian Levenberg-Marquardt SVD refinement on %d variables:'%(n)) 160 161 Lam = np.zeros((n,n)) 161 162 while icycle < maxcyc: … … 174 175 Amat /= Anorm 175 176 if Print: 176 print('initial chi^2 %.5g on %d obs.'%(chisq0,Nobs))177 G2fil.G2Print('initial chi^2 %.5g on %d obs.'%(chisq0,Nobs)) 177 178 chitol = ftol 178 179 while True: … … 183 184 except nl.LinAlgError: 184 185 psing = list(np.where(np.diag(nl.qr(Amatlam)[1]) < 1.e-14)[0]) 185 print ('ouch #1 bad SVD inversion; change parameterization for ',psing)186 G2fil.G2Print('ouch #1 bad SVD inversion; change parameterization for ',psing, mode='error') 186 187 return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-1}] 187 188 Xvec = np.inner(Ainv,Yvec) #solve … … 193 194 lam *= 10. 194 195 if Print: 195 print('new chi^2 %.5g on %d obs., %d SVD zeros ; matrix modification needed; lambda now %.1e' \196 G2fil.G2Print('new chi^2 %.5g on %d obs., %d SVD zeros ; matrix modification needed; lambda now %.1e' \ 196 197 %(chisq1,Nobs,Nzeros,lam)) 197 198 else: … … 200 201 break 201 202 if lam > 10.: 202 print ('ouch #3 chisq1 %g.4 stuck > chisq0 %g.4'%(chisq1,chisq0))203 G2fil.G2Print('ouch #3 chisq1 %g.4 stuck > chisq0 %g.4'%(chisq1,chisq0), mode='warn') 203 204 break 204 205 chitol *= 2 … … 206 207 deltaChi2 = (chisq0-chisq1)/chisq0 207 208 if Print: 208 print(' Cycle: %d, Time: %.2fs, Chi**2: %.5g for %d obs., Lambda: %.3g, Delta: %.3g'%(209 G2fil.G2Print(' Cycle: %d, Time: %.2fs, Chi**2: %.5g for %d obs., Lambda: %.3g, Delta: %.3g'%( 209 210 icycle,time.time()-time0,chisq1,Nobs,lamMax,deltaChi2)) 210 211 Histograms = args[0][0] … … 212 213 if deltaChi2 < ftol: 213 214 ifConverged = True 214 if Print: print("converged")215 if Print: G2fil.G2Print("converged") 215 216 break 216 217 icycle += 1 … … 224 225 try: 225 226 Bmat,Nzero = pinv(Amatlam,xtol) #Moore-Penrose inversion (via SVD) & count of zeros 226 if Print: print ('Found %d SVD zeros'%(Nzero))227 if Print: G2fil.G2Print('Found %d SVD zeros'%(Nzero), mode='warn') 227 228 Bmat = Bmat/Anorm 228 229 return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[],'SVD0':Nzero,'Converged': ifConverged, 'DelChi2':deltaChi2}] 229 230 except nl.LinAlgError: 230 print ('ouch #2 linear algebra error in making v-cov matrix')231 G2fil.G2Print('ouch #2 linear algebra error in making v-cov matrix', mode='error') 231 232 psing = [] 232 233 if maxcyc: … … 289 290 nfev = 0 290 291 if Print: 291 print(' Hessian SVD refinement on %d variables:'%(n))292 G2fil.G2Print(' Hessian SVD refinement on %d variables:'%(n)) 292 293 while icycle < maxcyc: 293 294 time0 = time.time() … … 304 305 Amat /= Anorm 305 306 if Print: 306 print('initial chi^2 %.5g'%(chisq0))307 G2fil.G2Print('initial chi^2 %.5g'%(chisq0)) 307 308 try: 308 309 Ainv,Nzeros = pinv(Amat,xtol) #do Moore-Penrose inversion (via SVD) 309 310 except nl.LinAlgError: 310 print ('ouch #1 bad SVD inversion; change parameterization')311 G2fil.G2Print('ouch #1 bad SVD inversion; change parameterization', mode='warn') 311 312 psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0]) 312 313 return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':0.,'psing':psing,'SVD0':-1}] … … 318 319 deltaChi2 = (chisq0-chisq1)/chisq0 319 320 if Print: 320 print(' Cycle: %d, Time: %.2fs, Chi**2: %.5g, Delta: %.3g'%(321 G2fil.G2Print(' Cycle: %d, Time: %.2fs, Chi**2: %.5g, Delta: %.3g'%( 321 322 icycle,time.time()-time0,chisq1,deltaChi2)) 322 323 Histograms = args[0][0] … … 324 325 if deltaChi2 < ftol: 325 326 ifConverged = True 326 if Print: print("converged")327 if Print: G2fil.G2Print("converged") 327 328 break 328 329 icycle += 1 … … 335 336 try: 336 337 Bmat,Nzero = pinv(Amat,xtol) #Moore-Penrose inversion (via SVD) & count of zeros 337 print ('Found %d SVD zeros'%(Nzero))338 G2fil.G2Print('Found %d SVD zeros'%(Nzero), mode='warn') 338 339 # Bmat = nl.inv(Amatlam); Nzeros = 0 339 340 Bmat = Bmat/Anorm … … 341 342 'SVD0':Nzero,'Converged': ifConverged, 'DelChi2':deltaChi2}] 342 343 except nl.LinAlgError: 343 print ('ouch #2 linear algebra error in making v-cov matrix')344 G2fil.G2Print('ouch #2 linear algebra error in making v-cov matrix', mode='error') 344 345 psing = [] 345 346 if maxcyc: … … 2741 2742 Boxes[box[0],box[1],box[2],Boxes[box[0],box[1],box[2],0]] = ib 2742 2743 except IndexError: 2743 print('too many atoms in box' )2744 G2fil.G2Print('Error: too many atoms in box' ) 2744 2745 continue 2745 2746 #Box content checks with errat.f $ erratv2.cpp ibox1 arrays … … 3010 3011 Values2Dict(parmDict, varyList, result[0]) 3011 3012 GOF = chisq/(len(result[2]['fvec'])-len(varyList)) #reduced chi^2 3012 print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],len(result[2]['fvec']),len(varyList)))3013 print ('refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc))3013 G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],len(result[2]['fvec']),len(varyList))) 3014 G2fil.G2Print ('refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)) 3014 3015 try: 3015 3016 sig = np.sqrt(np.diag(result[1])*GOF) 3016 3017 if np.any(np.isnan(sig)): 3017 print ('*** Least squares aborted - some invalid esds possible ***')3018 G2fil.G2Print ('*** Least squares aborted - some invalid esds possible ***', mode='error') 3018 3019 break #refinement succeeded - finish up! 3019 3020 except ValueError: #result[1] is None on singular matrix 3020 print ('**** Refinement failed - singular matrix ****')3021 G2fil.G2Print ('**** Refinement failed - singular matrix ****', mode='error') 3021 3022 return None 3022 3023 else: … … 3065 3066 generalData = data['General'] 3066 3067 if not generalData['Map']['MapType']: 3067 print ('**** ERROR - Fourier map not defined')3068 G2fil.G2Print ('**** ERROR - Fourier map not defined') 3068 3069 return 3069 3070 mapData = generalData['Map'] … … 3119 3120 mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho'])) 3120 3121 mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])] 3121 print ('Omit map time: %.4f no. elements: %d dimensions: %s'%(time.time()-time0,Fhkl.size,str(Fhkl.shape)))3122 G2fil.G2Print ('Omit map time: %.4f no. elements: %d dimensions: %s'%(time.time()-time0,Fhkl.size,str(Fhkl.shape))) 3122 3123 return mapData 3123 3124 … … 3185 3186 Fhkl[h,k,l] = complex(Fosq,0.) 3186 3187 rho = fft.fftn(fft.fftshift(Fhkl))/cell[6] 3187 print ('Fourier map time: %.4f no. elements: %d dimensions: %s'%(time.time()-time0,Fhkl.size,str(Fhkl.shape)))3188 G2fil.G2Print ('Fourier map time: %.4f no. elements: %d dimensions: %s'%(time.time()-time0,Fhkl.size,str(Fhkl.shape))) 3188 3189 mapData['Type'] = reflDict['Type'] 3189 3190 mapData['rho'] = np.real(rho) … … 3256 3257 mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho'])) 3257 3258 mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])] 3258 print ('Fourier map time: %.4f no. elements: %d dimensions: %s'%(time.time()-time0,Fhkl.size,str(Fhkl.shape)))3259 G2fil.G2Print ('Fourier map time: %.4f no. elements: %d dimensions: %s'%(time.time()-time0,Fhkl.size,str(Fhkl.shape))) 3259 3260 3260 3261 # map printing for testing purposes … … 3339 3340 i += 1 3340 3341 DH = np.array(DH) 3341 print (' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist)))3342 G2fil.G2Print (' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))) 3342 3343 Dphi = np.array(Dphi) 3343 3344 steps = np.array(hklShape) … … 3351 3352 chisq = np.min(Mmap) 3352 3353 DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape)) 3353 print (' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2]))3354 G2fil.G2Print (' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])) 3354 3355 # print (np.dot(DX,DH.T)+.5)%1.-Dphi 3355 3356 return DX … … 3444 3445 break 3445 3446 np.seterr(**old) 3446 print (' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size))3447 G2fil.G2Print (' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)) 3447 3448 CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))/10. #? to get on same scale as e-map 3448 print (' No.cycles = %d Residual Rcf =%8.3f%s Map size: %s'%(Ncyc,Rcf,'%',str(CErho.shape)))3449 G2fil.G2Print (' No.cycles = %d Residual Rcf =%8.3f%s Map size: %s'%(Ncyc,Rcf,'%',str(CErho.shape))) 3449 3450 roll = findOffset(SGData,A,CEhkl) #CEhkl needs to be just the observed set, not the full set! 3450 3451 … … 3505 3506 i += 1 3506 3507 DH = np.array(DH) 3507 print (' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist)))3508 G2fil.G2Print (' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))) 3508 3509 Dphi = np.array(Dphi) 3509 3510 steps = np.array(hklmShape) … … 3517 3518 chisq = np.min(Mmap) 3518 3519 DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape)) 3519 print (' map offset chi**2: %.3f, map offset: %d %d %d %d'%(chisq,DX[0],DX[1],DX[2],DX[3]))3520 G2fil.G2Print (' map offset chi**2: %.3f, map offset: %d %d %d %d'%(chisq,DX[0],DX[1],DX[2],DX[3])) 3520 3521 # print (np.dot(DX,DH.T)+.5)%1.-Dphi 3521 3522 return DX … … 3605 3606 break 3606 3607 np.seterr(**old) 3607 print (' Charge flip time: %.4f no. elements: %d'%(time.time()-time0,Ehkl.size))3608 G2fil.G2Print (' Charge flip time: %.4f no. elements: %d'%(time.time()-time0,Ehkl.size)) 3608 3609 CErho = np.real(fft.fftn(fft.fftshift(CEhkl[:,:,:,maxM+1])))/10. #? to get on same scale as e-map 3609 3610 SSrho = np.real(fft.fftn(fft.fftshift(CEhkl)))/10. #? ditto 3610 print (' No.cycles = %d Residual Rcf =%8.3f%s Map size: %s'%(Ncyc,Rcf,'%',str(CErho.shape)))3611 G2fil.G2Print (' No.cycles = %d Residual Rcf =%8.3f%s Map size: %s'%(Ncyc,Rcf,'%',str(CErho.shape))) 3611 3612 roll = findSSOffset(SGData,SSGData,A,CEhkl) #CEhkl needs to be just the observed set, not the full set! 3612 3613 … … 3778 3779 steps = np.array((3*[step,]),dtype='int32') 3779 3780 except KeyError: 3780 print ('**** ERROR - Fourier map not defined')3781 G2fil.G2Print ('**** ERROR - Fourier map not defined') 3781 3782 return peaks,mags 3782 3783 rhoMask = ma.array(rho,mask=(rho<contLevel)) … … 4436 4437 af = asarray(fqueue)*1.0 4437 4438 if retval == 5: 4438 print ('User terminated run; incomplete MC/SA')4439 G2fil.G2Print ('Error: User terminated run; incomplete MC/SA') 4439 4440 keepGoing = False 4440 4441 break … … 4443 4444 if abs(af[-1]-best_state.cost) > feps*10: 4444 4445 retval = 5 4445 print (" Warning: Cooled to %.4f > selected Tmin %.4f in %d steps"%(squeeze(last_state.cost),Tf,iters-1))4446 G2fil.G2Print (" Warning: Cooled to %.4f > selected Tmin %.4f in %d steps"%(squeeze(last_state.cost),Tf,iters-1)) 4446 4447 break 4447 4448 if (Tf is not None) and (schedule.T < Tf): … … 4453 4454 break 4454 4455 if (iters > maxiter): 4455 print (" Warning: Maximum number of iterations exceeded.")4456 G2fil.G2Print (" Warning: Maximum number of iterations exceeded.") 4456 4457 retval = 3 4457 4458 break … … 4473 4474 timelist.append(result[1]) 4474 4475 nsflist.append(result[2]) 4475 print (' MC/SA final fit: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1]))4476 G2fil.G2Print (' MC/SA final fit: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1])) 4476 4477 out_q.put(outlist) 4477 4478 out_t.put(timelist) … … 4879 4880 refs = np.array(refs).T 4880 4881 if start: 4881 print (' Minimum d-spacing used: %.2f No. reflections used: %d'%(MCSA['dmin'],nRef))4882 print (' Number of parameters varied: %d'%(len(varyList)))4882 G2fil.G2Print (' Minimum d-spacing used: %.2f No. reflections used: %d'%(MCSA['dmin'],nRef)) 4883 G2fil.G2Print (' Number of parameters varied: %d'%(len(varyList))) 4883 4884 start = False 4884 4885 parmDict['sumFosq'] = sumFosq … … 4902 4903 lower=lower, upper=upper, slope=MCSA['log slope'],ranStart=MCSA.get('ranStart',False), 4903 4904 ranRange=MCSA.get('ranRange',10.)/100.,autoRan=MCSA.get('autoRan',False),dlg=pgbar) 4904 print (' Acceptance rate: %.2f%% MCSA residual: %.2f%%'%(100.*results[5]/results[3],100.*results[1]))4905 G2fil.G2Print (' Acceptance rate: %.2f%% MCSA residual: %.2f%%'%(100.*results[5]/results[3],100.*results[1])) 4905 4906 results = so.minimize(mcsaCalc,results[0],method='L-BFGS-B',args=(refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict), 4906 4907 bounds=bounds,)
Note: See TracChangeset
for help on using the changeset viewer.