Changeset 4789 for trunk/GSASIImath.py
- Timestamp:
- Feb 1, 2021 1:36:50 PM (2 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIImath.py
r4770 r4789 126 126 127 127 def setHcorr(info,Amat,xtol,problem=False): 128 '''Find & report high correlation terms in covariance matrix 129 ''' 128 130 Adiag = np.sqrt(np.diag(Amat)) 129 131 if np.any(np.abs(Adiag) < 1.e-14): raise G2NormException # test for any hard singularities … … 141 143 info['Hcorr'] = [(i,j,AcovUp[i,j]) for i,j in zip(*np.where(AcovUp > m))] 142 144 return Bmat,Nzeros 143 145 146 def setSVDwarn(info,Amat,Nzeros,indices): 147 '''Find & report terms causing SVN zeros 148 ''' 149 if Nzeros == 0: return 150 d = np.abs(np.diag(nl.qr(Amat)[1])) 151 svdsing = list(np.where(d < 1.e-14)[0]) 152 if len(svdsing) < Nzeros: # try to get the Nzeros worst terms 153 svdsing = list(np.where(d <= sorted(d)[Nzeros-1])[0]) 154 155 156 if not len(svdsing): # make sure at least the worst term is shown 157 svdsing = [np.argmin(d)] 158 info['SVDsing'] = [indices[i] for i in svdsing] 159 144 160 def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.e-6, maxcyc=0,lamda=-3,Print=False,refPlotUpdate=None): 145 161 ''' … … 187 203 ifConverged = False 188 204 deltaChi2 = -10. 189 x0 = np.array(x0, ndmin=1, dtype=np.float64) #might be redundant? 205 x0 = np.array(x0, ndmin=1, dtype=np.float64) # make sure that x0 float 1-D 206 # array (in case any parameters were set to int) 190 207 n = len(x0) 191 208 if type(args) != type(()): … … 193 210 194 211 icycle = 0 212 AmatAll = None # changed only if cycles > 0 195 213 lamMax = lam = 0 # start w/o Marquardt and add it in if needed 196 214 nfev = 0 … … 221 239 nfev += 1 222 240 chisq0 = np.sum(M**2) 223 Yvec,AmatAll = Hess(x0,*args) # compute hessian & vectors with all variables 241 YvecAll,AmatAll = Hess(x0,*args) # compute hessian & vectors with all variables 242 Yvec = copy.copy(YvecAll) 224 243 Amat = copy.copy(AmatAll) 225 244 Xvec = copy.copy(XvecAll) … … 257 276 info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros} 258 277 info['psing'] = [i for i in range(n) if i not in indices] 259 260 278 info['msg'] = 'SVD inversion failure\n' 261 279 return [x0,None,info] 262 if Nzeros: # remove bad vars from the matrix263 d = np.abs(np.diag(nl.qr(Amatlam)[1]))264 psing = list(np.where(d < 1.e-14)[0])265 if not len(psing): # make sure at least the worst term is removed266 psing = [np.argmin(d)]267 G2fil.G2Print('{} SVD Zeros: dropping terms for for variable(s) #{}'.268 format(Nzeros,psing), mode='warn')269 Amat, indices, Xvec, Yvec, Adiag = dropTerms(psing,Amat, indices, Xvec, Yvec, Adiag)270 Amatlam = Amat*(1.+np.eye(Amat.shape[0])*lam)271 Ainv,nz = pinv(Amatlam,xtol) #do Moore-Penrose inversion (via SVD)272 if nz > 0: G2fil.G2Print('Note: there are {} new SVD Zeros after drop'.format(nz),273 mode='warn')274 280 Xvec = np.inner(Ainv,Yvec)/Adiag #solve for LS terms 275 281 XvecAll[indices] = Xvec # expand … … 295 301 G2fil.G2Print(Msg.msg,mode='warn') 296 302 info['msg'] = Msg.msg + '\n\n' 303 setSVDwarn(info,Amatlam,Nzeros,indices) 297 304 return [x0,None,info] 298 305 nfev += 1 … … 311 318 try: # report highly correlated parameters from full Hessian, if we can 312 319 info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax, 313 'Converged':False, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00} 320 'Converged':False, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 321 'chisq0':chisq00, 'Ouch#4':True} 314 322 info['psing'] = [i for i in range(n) if i not in indices] 315 323 Bmat,Nzeros = setHcorr(info,AmatAll,xtol,problem=True) 316 324 info['SVD0'] = Nzeros 325 setSVDwarn(info,Amatlam,Nzeros,indices) 317 326 return [x0,Bmat,info] 318 327 except Exception as Msg: … … 331 340 if Print and n-len(indices): 332 341 G2fil.G2Print( 333 'Cycle %d: %.2fs Chi2: %.5g; Obs: %d; Lam: %.3g Del: %.3g; drop=%d '%334 (icycle,time.time()-time0,chisq1,Nobs,lamMax,deltaChi2,n-len(indices) ))342 'Cycle %d: %.2fs Chi2: %.5g; Obs: %d; Lam: %.3g Del: %.3g; drop=%d, SVD=%d'% 343 (icycle,time.time()-time0,chisq1,Nobs,lamMax,deltaChi2,n-len(indices),Nzeros)) 335 344 elif Print: 336 345 G2fil.G2Print( 337 'Cycle %d: %.2fs, Chi**2: %.5g for %d obs., Lambda: %.3g, Delta: %.3g '%338 (icycle,time.time()-time0,chisq1,Nobs,lamMax,deltaChi2 ))346 'Cycle %d: %.2fs, Chi**2: %.5g for %d obs., Lambda: %.3g, Delta: %.3g, SVD=%d'% 347 (icycle,time.time()-time0,chisq1,Nobs,lamMax,deltaChi2,Nzeros)) 339 348 Histograms = args[0][0] 340 349 if refPlotUpdate is not None: refPlotUpdate(Histograms,icycle) # update plot … … 355 364 info = {'num cyc':icycle,'fvec':M2,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-2} 356 365 info['msg'] = Msg.msg + '\n' 366 setSVDwarn(info,Amatlam,Nzeros,indices) 357 367 return [x0,None,info] 358 368 chisqf = np.sum(M**2) # ending chi**2 359 Yvec,Amat = Hess(x0,*args) # we could save some time and use the last Hessian from the last refinement cycle360 369 psing_prev = [i for i in range(n) if i not in indices] # save dropped vars 370 if AmatAll is None: # Save some time and use Hessian from the last refinement cycle 371 Yvec,Amat = Hess(x0,*args) 372 else: 373 Yvec = copy.copy(YvecAll) 374 Amat = copy.copy(AmatAll) 361 375 indices = range(n) 362 376 info = {} … … 366 380 'Converged':ifConverged, 'DelChi2':deltaChi2, 'chisq0':chisq00}) 367 381 if icycle > 0: info.update({'Xvec':XvecAll}) 382 setSVDwarn(info,Amatlam,Nzeros,indices) 383 # expand Bmat by filling with zeros if columns have been dropped 384 if len(psing_prev): 385 ins = [j-i for i,j in enumerate(psing_prev)] 386 Bmat = np.insert(np.insert(Bmat,ins,0,1),ins,0,0) 368 387 return [x0,Bmat,info] 369 388 except nl.LinAlgError: … … 392 411 info['psing'] = [i for i in range(n) if i not in indices] 393 412 return [x0,None,info] 394 if Nzeros: # sometimes SVD zeros show up when lam=0395 d = np.abs(np.diag(nl.qr(Amat)[1]))396 psing = list(np.where(d < 1.e-14)[0])397 if len(psing) < Nzeros:398 psing = list(np.where(d <= sorted(d)[Nzeros-1])[0])399 if Print:400 G2fil.G2Print('Found %d SVD zeros w/o Lambda. '+401 'Likely problem with variable(s) #%s'%(Nzeros,psing), mode='warn')402 Amat, indices, Yvec = dropTerms(psing, Amat, indices, Yvec)403 413 # expand Bmat by filling with zeros if columns have been dropped 404 414 psing = [i for i in range(n) if i not in indices] … … 409 419 'Converged':ifConverged, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00}) 410 420 info['psing'] = [i for i in range(n) if i not in indices] 421 setSVDwarn(info,Amat,Nzeros,indices) 411 422 return [x0,Bmat,info] 412 423
Note: See TracChangeset
for help on using the changeset viewer.