Changeset 4690 for trunk/GSASIImath_new.py
 Timestamp:
 Dec 30, 2020 5:52:57 PM (2 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIImath_new.py
r4683 r4690 56 56 ##### Hessian leastsquares LevenbergMarquardt routine 57 57 ################################################################################ 58 58 class G2NormException(Exception): pass 59 59 60 def pinv(a, rcond=1e15 ): 60 61 ''' … … 124 125 np.delete(yvector,bad), 125 126 np.delete(indices,bad)) 127 128 def setHcorr(info,Amat,xtol,problem=False): 129 Adiag = np.sqrt(np.diag(Amat)) 130 if np.any(np.abs(Adiag) < 1.e14): raise G2NormException # test for any hard singularities 131 Anorm = np.outer(Adiag,Adiag) 132 Amat = Amat/Anorm 133 Bmat,Nzeros = pinv(Amat,xtol) #MoorePenrose inversion (via SVD) & count of zeros 134 Bmat = Bmat/Anorm 135 sig = np.sqrt(np.diag(Bmat)) 136 xvar = np.outer(sig,np.ones_like(sig)) 137 AcovUp = abs(np.triu(np.divide(np.divide(Bmat,xvar),xvar.T),1)) # elements above diagonal 138 if Nzeros or problem: # something is wrong, so report what is found 139 m = min(0.99,0.99*np.amax(AcovUp)) 140 else: 141 m = max(0.95,0.99*np.amax(AcovUp)) 142 info['Hcorr'] = [(i,j,AcovUp[i,j]) for i,j in zip(*np.where(AcovUp > m))] 143 return Bmat 126 144 127 145 def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e8,xtol=1.e6, maxcyc=0,lamda=3,Print=False,refPlotUpdate=None): … … 195 213 if Print: 196 214 G2fil.G2Print('initial chi^2 %.5g with %d obs.'%(chisq00,Nobs)) 215 if n == 0: 216 info = {'num cyc':0,'fvec':M2,'nfev':1,'lamMax':0,'psing':[],'SVD0':0} 217 info['msg'] = 'no variables, no refinement\n' 218 return [x0,None,info] 197 219 while icycle < maxcyc: 198 220 M = M2 … … 208 230 Adiag = np.sqrt(np.diag(Amat)) 209 231 psing = np.where(np.abs(Adiag) < 1.e14)[0] # remove any hard singularities 210 if psing:232 if len(psing): 211 233 G2fil.G2Print('ouch #1 dropping singularities for variable(s) #{}'.format(psing), mode='warn') 212 234 Amat, Xvec, Yvec, indices = dropTerms(Amat, Xvec, Yvec, indices, psing) 213 235 Adiag = np.sqrt(np.diag(Amat)) 214 Anorm = np.outer(Adiag [indices],Adiag[indices]) # normalize matrix & vector215 Yvec /= Adiag [indices]236 Anorm = np.outer(Adiag,Adiag) # normalize matrix & vector 237 Yvec /= Adiag 216 238 Amat /= Anorm 217 239 chitol = ftol … … 226 248 d = np.abs(np.diag(nl.qr(Amatlam)[1])) 227 249 psing = list(np.where(d < 1.e14)[0]) 228 if not psing: # make sure at least the worst term is removed250 if not len(psing): # make sure at least the worst term is removed 229 251 psing = [np.argmin(d)] 230 252 G2fil.G2Print('ouch #2 bad SVD inversion; dropping terms for for variable(s) #{}'. 231 253 format(psing), mode='warn') 232 254 Amatlam, Xvec, Yvec, indices = dropTerms(Amatlam, Xvec, Yvec, indices, psing) 233 Anorm = np.outer(Adiag[indices],Adiag[indices])234 255 if loops < maxdrop: continue 235 256 G2fil.G2Print('giving up with ouch #2', mode='error') … … 240 261 d = np.abs(np.diag(nl.qr(Amatlam)[1])) 241 262 psing = list(np.where(d < 1.e14)[0]) 242 if not psing: # make sure at least the worst term is removed263 if not len(psing): # make sure at least the worst term is removed 243 264 psing = [np.argmin(d)] 244 265 G2fil.G2Print('{} SVD Zeros: dropping terms for for variable(s) #{}'. … … 246 267 Amat, Xvec, Yvec, indices = dropTerms(Amat, Xvec, Yvec, indices, psing) 247 268 Amatlam = Amat*(1.+np.eye(Amat.shape[0])*lam) 248 Anorm = np.outer(Adiag[indices],Adiag[indices])249 269 Ainv,Nzeros = pinv(Amatlam,xtol) #do MoorePenrose inversion (via SVD) 250 Xvec = np.inner(Ainv,Yvec)/Adiag[indices] #solve for LS terms 270 Adiag = np.sqrt(np.diag(Amat)) 271 Xvec = np.inner(Ainv,Yvec)/Adiag #solve for LS terms 251 272 XvecAll[indices] = Xvec # expand 252 273 try: … … 261 282 G2fil.G2Print(Msg.msg,mode='warn') 262 283 # dropping variables does not seem to help 263 if psing:284 if len(psing): 264 285 G2fil.G2Print('ouch #3 unable to evaluate objective function;'+ 265 286 '\n\tdropping terms for variable(s) #{}'.format(psing), … … 274 295 info['psing'] = [i for i in range(n) if i not in indices] 275 296 try: # report highly correlated parameters from full Hessian, if we can 276 Adiag = np.sqrt(np.diag(AmatAll)) 277 Anorm = np.outer(Adiag,Adiag) 278 if np.where(np.abs(Adiag) < 1.e14)[0]: raise Exception # test for any hard singularities 279 Amat = AmatAll/Anorm 280 Bmat,Nzeros = pinv(Amat,xtol) #MoorePenrose inversion (via SVD) & count of zeros 281 Bmat = Bmat/Anorm 282 sig = np.sqrt(np.diag(Bmat)) 283 xvar = np.outer(sig,np.ones_like(sig)) 284 AcovUp = abs(np.triu(np.divide(np.divide(Bmat,xvar),xvar.T),1)) # elements above diagonal 285 m = min(0.99,0.99*np.amax(AcovUp)) 286 info['Hcorr'] = [(i,j,AcovUp[i,j]) for i,j in zip(*np.where(AcovUp > m))] 287 except: 288 pass 297 setHcorr(info,AmatAll,xtol,problem=True) 298 except nl.LinAlgError: 299 G2fil.G2Print('Warning: Hessian too illconditioned to get full covariance matrix') 300 except G2NormException: 301 G2fil.G2Print('Warning: Hessian normalization problem') 302 except Exception as Msg: 303 if not hasattr(Msg,'msg'): Msg.msg = str(Msg) 304 G2fil.G2Print('Error determining highly correlated vars',mode='warn') 305 G2fil.G2Print(Msg.msg,mode='warn') 289 306 info['msg'] = Msg.msg + '\n' 290 307 return [x0,None,info] … … 298 315 '\tincreasing Marquardt lambda to %.1e')%(chisq1,Nobs,Nzeros,lam)) 299 316 if lam > 10.: 300 G2fil.G2Print('ouch #4 stuck: chisqnew % g.4 > chisq0 %g.4 with lambda %g.1'%317 G2fil.G2Print('ouch #4 stuck: chisqnew %.4g > chisq0 %.4g with lambda %.1g'% 301 318 (chisq1,chisq0,lam), mode='warn') 319 try: # report highly correlated parameters from full Hessian, if we can 320 info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros, 321 'Converged':False, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00} 322 info['psing'] = [i for i in range(n) if i not in indices] 323 Bmat = setHcorr(info,AmatAll,xtol,problem=True) 324 return [x0,Bmat,info] 325 except Exception as Msg: 326 if not hasattr(Msg,'msg'): Msg.msg = str(Msg) 327 G2fil.G2Print('Error determining highly correlated vars',mode='warn') 328 G2fil.G2Print(Msg.msg,mode='warn') 302 329 maxcyc = 1 # no more cycles 303 330 break 304 331 chitol *= 2 305 332 else: # refinement succeeded 306 x0 [indices] += Xvec333 x0 += XvecAll 307 334 lam /= 10. # drop lam on next cycle 308 335 break … … 343 370 info = {} 344 371 try: # report highly correlated parameters from full Hessian, if we can 345 Adiag = np.sqrt(np.diag(Amat)) 346 Anorm = np.outer(Adiag,Adiag) 347 if np.where(np.abs(Adiag) < 1.e14)[0]: raise Exception # test for any hard singularities 348 AmatN = Amat/Anorm 349 Bmat,Nzeros = pinv(AmatN,xtol) #MoorePenrose inversion (via SVD) & count of zeros 350 Bmat = Bmat/Anorm 351 sig = np.sqrt(np.diag(Bmat)) 352 xvar = np.outer(sig,np.ones_like(sig)) 353 AcovUp = abs(np.triu(np.divide(np.divide(Bmat,xvar),xvar.T),1)) # elements above diagonal 354 m = min(0.99,0.99*np.amax(AcovUp)) 355 info['Hcorr'] = [(i,j,AcovUp[i,j]) for i,j in zip(*np.where(AcovUp > m))] 372 Bmat = setHcorr(info,Amat,xtol,problem=False) 356 373 info.update({'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros,'psing':psing_prev, 357 374 'Converged':ifConverged, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00}) … … 359 376 except nl.LinAlgError: 360 377 G2fil.G2Print('Warning: Hessian too illconditioned to get full covariance matrix') 378 except G2NormException: 379 G2fil.G2Print('Warning: Hessian normalization problem') 361 380 except Exception as Msg: 362 381 if not hasattr(Msg,'msg'): Msg.msg = str(Msg) 382 G2fil.G2Print('Error determining highly correlated vars',mode='warn') 363 383 G2fil.G2Print(Msg.msg,mode='warn') 364 G2fil.G2Print('Warning: Hessian normalization problem') 365 # matrix inversion failed drop previously removed variables 366 Xvec = np.zeros(n) # dummy variable 367 Amat, Xvec, Yvec, indices = dropTerms(Amat, Xvec, Yvec, indices, psing_prev) 384 # matrix above inversion failed, drop previously removed variables 385 dummy = np.zeros(n) 386 Amat, Xvec, Yvec, indices = dropTerms(Amat, dummy, Yvec, indices, psing_prev) 368 387 Adiag = np.sqrt(np.diag(Amat)) 369 388 Anorm = np.outer(Adiag,Adiag) … … 375 394 G2fil.G2Print('ouch #5 linear algebra error in making final vcov matrix', mode='error') 376 395 psing = list(np.where(np.abs(np.diag(nl.qr(Amat)[1])) < 1.e14)[0]) 377 if not psing: # make sure at least the worst is returned396 if not len(psing): # make sure at least the worst is returned 378 397 psing = [np.argmin(d)] 379 398 Amat, Xvec, Yvec, indices = dropTerms(Amat, Xvec, Yvec, indices, psing) … … 391 410 Amat, Xvec, Yvec, indices = dropTerms(Amat, Xvec, Yvec, indices, psing) 392 411 # expand Bmat by filling with zeros 393 if psing:412 if len(psing): 394 413 Bmat = np.insert(np.insert(Bmat,psing,0,1),psing,0,0) 395 414 info.update({'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros,
Note: See TracChangeset
for help on using the changeset viewer.