Changeset 4709 for trunk/GSASIImath.py
 Timestamp:
 Jan 4, 2021 10:41:46 AM (2 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIImath.py
r4688 r4709 56 56 ##### Hessian leastsquares LevenbergMarquardt routine 57 57 ################################################################################ 58 58 class G2NormException(Exception): pass 59 59 60 def pinv(a, rcond=1e15 ): 60 61 ''' … … 101 102 s = np.where(s>cutoff,1./s,0.) 102 103 nzero = s.shape[0]np.count_nonzero(s) 103 # res = np.dot(np.transpose(vt), np.multiply(s[:, np.newaxis], np.transpose(u)))104 104 res = np.dot(vt.T,s[:,nxs]*u.T) 105 105 return res,nzero 106 106 107 def dropTerms(bad, hessian, indices, *vectors): 108 '''Remove the 'bad' terms from the Hessian and vector 109 110 :param tuple bad: a list of variable (row/column) numbers that should be 111 removed from the hessian and vector. Example: (0,3) removes the 1st and 112 4th column/row 113 :param np.array hessian: a square matrix of length n x n 114 :param np.array indices: the indices of the leastsquares vector of length n 115 referenced to the initial variable list; as this routine is called 116 multiple times, more terms may be removed from this list 117 :param additionalargs: various leastsquares model values, length n 118 119 :returns: hessian, indices, vector0, vector1,... where the lengths are 120 now n' x n' and n', with n' = n  len(bad) 121 ''' 122 out = [np.delete(np.delete(hessian,bad,1),bad,0),np.delete(indices,bad)] 123 for v in vectors: 124 out.append(np.delete(v,bad)) 125 return out 126 127 def setHcorr(info,Amat,xtol,problem=False): 128 Adiag = np.sqrt(np.diag(Amat)) 129 if np.any(np.abs(Adiag) < 1.e14): raise G2NormException # test for any hard singularities 130 Anorm = np.outer(Adiag,Adiag) 131 Amat = Amat/Anorm 132 Bmat,Nzeros = pinv(Amat,xtol) #MoorePenrose inversion (via SVD) & count of zeros 133 Bmat = Bmat/Anorm 134 sig = np.sqrt(np.diag(Bmat)) 135 xvar = np.outer(sig,np.ones_like(sig)) 136 AcovUp = abs(np.triu(np.divide(np.divide(Bmat,xvar),xvar.T),1)) # elements above diagonal 137 if Nzeros or problem: # something is wrong, so report what is found 138 m = min(0.99,0.99*np.amax(AcovUp)) 139 else: 140 m = max(0.95,0.99*np.amax(AcovUp)) 141 info['Hcorr'] = [(i,j,AcovUp[i,j]) for i,j in zip(*np.where(AcovUp > m))] 142 return Bmat,Nzeros 143 107 144 def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e8,xtol=1.e6, maxcyc=0,lamda=3,Print=False,refPlotUpdate=None): 108 145 ''' … … 138 175 residual standard deviation to get the covariance of the 139 176 parameter estimates  see curve_fit. 140 * infodict : dict 141 a dictionary of optional outputs with the keys: 177 * infodict : dict, a dictionary of optional outputs with the keys: 142 178 143 179 * 'fvec' : the function evaluated at the output 144 180 * 'num cyc': 145 * 'nfev': 181 * 'nfev': number of objective function evaluation calls 146 182 * 'lamMax': 147 * 'psing': 148 * 'SVD0': 149 150 ''' 151 183 * 'psing': list of variable variables that have been removed from the refinement 184 * 'SVD0': 1 for singlar matrix, 2 for objective function exception, Nzeroes = # of SVD 0's 185 * 'Hcorr': list entries (i,j,c) where i & j are of highly correlated variables & c is correlation coeff. 186 ''' 152 187 ifConverged = False 153 188 deltaChi2 = 10. … … 158 193 159 194 icycle = 0 160 One = np.ones((n,n)) 161 lam = 10.**lamda 162 lamMax = lam 195 lamMax = lam = 0 # start w/o Marquardt and add it in if needed 163 196 nfev = 0 164 197 if Print: 165 198 G2fil.G2Print(' Hessian LevenbergMarquardt SVD refinement on %d variables:'%(n)) 166 Lam = np.zeros((n,n)) 167 Xvec = np.zeros(len(x0)) 168 chisq00 = None 199 XvecAll = np.zeros(n) 200 try: 201 M2 = func(x0,*args) 202 except Exception as Msg: 203 if not hasattr(Msg,'msg'): Msg.msg = str(Msg) 204 G2fil.G2Print('\nouch #0 unable to evaluate initial objective function\nCheck for an invalid parameter value',mode='error') 205 G2fil.G2Print('Use Calculate/"View LS parms" to list all parameter values\n',mode='warn') 206 G2fil.G2Print('Error message: '+Msg.msg,mode='warn') 207 raise Exception('HessianLSQ  ouch #0: look for invalid parameter (see console)') 208 chisq00 = np.sum(M2**2) # starting Chi**2 209 Nobs = len(M2) 210 if Print: 211 G2fil.G2Print('initial chi^2 %.5g with %d obs.'%(chisq00,Nobs)) 212 if n == 0: 213 info = {'num cyc':0,'fvec':M2,'nfev':1,'lamMax':0,'psing':[],'SVD0':0} 214 info['msg'] = 'no variables: skipping refinement\n' 215 return [x0,None,info] 169 216 while icycle < maxcyc: 217 M = M2 170 218 time0 = time.time() 171 M = func(x0,*args)172 Nobs = len(M)173 219 nfev += 1 174 220 chisq0 = np.sum(M**2) 175 if chisq00 is None: chisq00 = chisq0 176 Yvec,Amat = Hess(x0,*args) 221 Yvec,AmatAll = Hess(x0,*args) # compute hessian & vectors with all variables 222 Amat = copy.copy(AmatAll) 223 Xvec = copy.copy(XvecAll) 224 # we could remove vars that were dropped in previous cycles here (use indices), but for 225 # now let's reset them each cycle in case the singularities go away 226 indices = range(n) 177 227 Adiag = np.sqrt(np.diag(Amat)) 178 #psing = np.where(np.abs(Adiag) < 1.e14,True,False) 179 psing = np.where(np.abs(Adiag) < 1.e14)[0] 180 if np.any(psing): #hard singularity in matrix 181 return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':1}] 182 Anorm = np.outer(Adiag,Adiag) 228 psing = np.where(np.abs(Adiag) < 1.e14)[0] # find any hard singularities 229 if len(psing): 230 G2fil.G2Print('ouch #1 dropping singularities for variable(s) #{}'.format( 231 psing), mode='warn') 232 Amat, indices, Xvec, Yvec, Adiag = dropTerms(psing,Amat, indices, Xvec, Yvec, Adiag) 233 Anorm = np.outer(Adiag,Adiag) # normalize matrix & vector 183 234 Yvec /= Adiag 184 235 Amat /= Anorm 185 if Print:186 G2fil.G2Print('initial chi^2 %.5g on %d obs.'%(chisq0,Nobs))187 236 chitol = ftol 188 while True: 189 Lam = np.eye(Amat.shape[0])*lam 190 Amatlam = Amat*(One+Lam) 237 maxdrop = 3 # max number of drop var attempts 238 loops = 0 239 while True: # loop as we increase lamda and possibly drop terms 240 lamMax = max(lamMax,lam) 241 Amatlam = Amat*(1.+np.eye(Amat.shape[0])*lam) 191 242 try: 192 Ainv,Nzeros = pinv(Amatlam,xtol) # do MoorePenrose inversion (via SVD)243 Ainv,Nzeros = pinv(Amatlam,xtol) # MoorePenrose SVD inversion 193 244 except nl.LinAlgError: 194 psing = list(np.where(np.abs(np.diag(nl.qr(Amatlam)[1])) < 1.e14)[0]) 195 G2fil.G2Print('ouch #1 bad SVD inversion; change parameterization for ',psing, mode='error') 196 return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':1}] 197 Xvec = np.inner(Ainv,Yvec) #solve 198 Xvec /= Adiag 199 M2 = func(x0+Xvec,*args) 245 loops += 1 246 d = np.abs(np.diag(nl.qr(Amatlam)[1])) 247 psing = list(np.where(d < 1.e14)[0]) 248 if not len(psing): # make sure at least the worst term is removed 249 psing = [np.argmin(d)] 250 G2fil.G2Print('ouch #2 bad SVD inversion; dropping terms for for variable(s) #{}'. 251 format(psing), mode='warn') 252 Amat, indices, Xvec, Yvec, Adiag = dropTerms(psing,Amat, indices, Xvec, Yvec, Adiag) 253 if loops < maxdrop: continue # try again, same lam but fewer vars 254 G2fil.G2Print('giving up with ouch #2', mode='error') 255 info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros} 256 info['psing'] = [i for i in range(n) if i not in indices] 257 258 info['msg'] = 'SVD inversion failure\n' 259 return [x0,None,info] 260 if Nzeros: # remove bad vars from the matrix 261 d = np.abs(np.diag(nl.qr(Amatlam)[1])) 262 psing = list(np.where(d < 1.e14)[0]) 263 if not len(psing): # make sure at least the worst term is removed 264 psing = [np.argmin(d)] 265 G2fil.G2Print('{} SVD Zeros: dropping terms for for variable(s) #{}'. 266 format(Nzeros,psing), mode='warn') 267 Amat, indices, Xvec, Yvec, Adiag = dropTerms(psing,Amat, indices, Xvec, Yvec, Adiag) 268 Amatlam = Amat*(1.+np.eye(Amat.shape[0])*lam) 269 Ainv,nz = pinv(Amatlam,xtol) #do MoorePenrose inversion (via SVD) 270 if nz > 0: G2fil.G2Print('Note: there are {} new SVD Zeros after drop'.format(nz), 271 mode='warn') 272 Xvec = np.inner(Ainv,Yvec)/Adiag #solve for LS terms 273 XvecAll[indices] = Xvec # expand 274 try: 275 M2 = func(x0+XvecAll,*args) 276 except Exception as Msg: 277 if not hasattr(Msg,'msg'): Msg.msg = str(Msg) 278 G2fil.G2Print(Msg.msg,mode='warn') 279 loops += 1 280 d = np.abs(np.diag(nl.qr(Amatlam)[1])) 281 G2fil.G2Print('ouch #3 unable to evaluate objective function;',mode='error') 282 info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros} 283 info['psing'] = [i for i in range(n) if i not in indices] 284 try: # try to report highly correlated parameters from full Hessian 285 setHcorr(info,AmatAll,xtol,problem=True) 286 except nl.LinAlgError: 287 G2fil.G2Print('Warning: Hessian too illconditioned to get full covariance matrix') 288 except G2NormException: 289 G2fil.G2Print('Warning: Hessian normalization problem') 290 except Exception as Msg: 291 if not hasattr(Msg,'msg'): Msg.msg = str(Msg) 292 G2fil.G2Print('Error determining highly correlated vars',mode='warn') 293 G2fil.G2Print(Msg.msg,mode='warn') 294 info['msg'] = Msg.msg + '\n\n' 295 return [x0,None,info] 200 296 nfev += 1 201 297 chisq1 = np.sum(M2**2) 202 298 if chisq1 > chisq0*(1.+chitol): #TODO put Alan Coehlo's criteria for lambda here? 203 lam *= 10. 299 if lam == 0: 300 lam = 10.**lamda # set to initial Marquardt term 301 else: 302 lam *= 10. 204 303 if Print: 205 G2fil.G2Print('new chi^2 %.5g on %d obs., %d SVD zeros ; matrix modification needed; lambda now %.1e' \ 206 %(chisq1,Nobs,Nzeros,lam)) 207 else: 208 x0 += Xvec 209 lam /= 10. 304 G2fil.G2Print(('divergence: chi^2 %.5g on %d obs. (%d SVD zeros)\n'+ 305 '\tincreasing Marquardt lambda to %.1e')%(chisq1,Nobs,Nzeros,lam)) 306 if lam > 10.: 307 G2fil.G2Print('ouch #4 stuck: chisqnew %.4g > chisq0 %.4g with lambda %.1g'% 308 (chisq1,chisq0,lam), mode='warn') 309 try: # report highly correlated parameters from full Hessian, if we can 310 info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax, 311 'Converged':False, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00} 312 info['psing'] = [i for i in range(n) if i not in indices] 313 Bmat,Nzeros = setHcorr(info,AmatAll,xtol,problem=True) 314 info['SVD0'] = Nzeros 315 return [x0,Bmat,info] 316 except Exception as Msg: 317 if not hasattr(Msg,'msg'): Msg.msg = str(Msg) 318 G2fil.G2Print('Error determining highly correlated vars',mode='warn') 319 G2fil.G2Print(Msg.msg,mode='warn') 320 maxcyc = 1 # no more cycles 321 break 322 chitol *= 2 323 else: # refinement succeeded 324 x0 += XvecAll 325 lam /= 10. # drop lam on next cycle 210 326 break 211 if lam > 10.: 212 G2fil.G2Print('ouch #3 chisq1 %g.4 stuck > chisq0 %g.4'%(chisq1,chisq0), mode='warn') 213 break 214 chitol *= 2 215 lamMax = max(lamMax,lam) 327 # complete current cycle 216 328 deltaChi2 = (chisq0chisq1)/chisq0 217 if Print: 218 G2fil.G2Print(' Cycle: %d, Time: %.2fs, Chi**2: %.5g for %d obs., Lambda: %.3g, Delta: %.3g'%( 219 icycle,time.time()time0,chisq1,Nobs,lamMax,deltaChi2)) 329 if Print and nlen(indices): 330 G2fil.G2Print( 331 'Cycle %d: %.2fs Chi2: %.5g; Obs: %d; Lam: %.3g Del: %.3g; drop=%d'% 332 (icycle,time.time()time0,chisq1,Nobs,lamMax,deltaChi2,nlen(indices))) 333 elif Print: 334 G2fil.G2Print( 335 'Cycle %d: %.2fs, Chi**2: %.5g for %d obs., Lambda: %.3g, Delta: %.3g'% 336 (icycle,time.time()time0,chisq1,Nobs,lamMax,deltaChi2)) 220 337 Histograms = args[0][0] 221 338 if refPlotUpdate is not None: refPlotUpdate(Histograms,icycle) # update plot … … 225 342 break 226 343 icycle += 1 227 M = func(x0,*args)344 # refinement complete, compute Covariance matrix w/o LevenbergMarquardt 228 345 nfev += 1 229 Yvec,Amat = Hess(x0,*args) 346 try: 347 M = func(x0,*args) 348 except Exception as Msg: 349 if not hasattr(Msg,'msg'): Msg.msg = str(Msg) 350 G2fil.G2Print(Msg.msg,mode='warn') 351 G2fil.G2Print('ouch #5 final objective function reeval failed',mode='error') 352 psing = [i for i in range(n) if i not in indices] 353 info = {'num cyc':icycle,'fvec':M2,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':2} 354 info['msg'] = Msg.msg + '\n' 355 return [x0,None,info] 356 chisqf = np.sum(M**2) # ending chi**2 357 Yvec,Amat = Hess(x0,*args) # we could save some time and use the last Hessian from the last refinement cycle 358 psing_prev = [i for i in range(n) if i not in indices] # save dropped vars 359 indices = range(n) 360 info = {} 361 try: # report highly correlated parameters from full Hessian, if we can 362 Bmat,Nzeros = setHcorr(info,Amat,xtol,problem=False) 363 info.update({'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros,'psing':psing_prev, 364 'Converged':ifConverged, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00}) 365 return [x0,Bmat,info] 366 except nl.LinAlgError: 367 G2fil.G2Print('Warning: Hessian too illconditioned to get full covariance matrix') 368 except G2NormException: 369 G2fil.G2Print('Warning: Hessian normalization problem') 370 except Exception as Msg: 371 if not hasattr(Msg,'msg'): Msg.msg = str(Msg) 372 G2fil.G2Print('Error determining highly correlated vars',mode='warn') 373 G2fil.G2Print(Msg.msg,mode='warn') 374 # matrix above inversion failed, drop previously removed variables & try again 375 Amat, indices, Yvec = dropTerms(psing_prev, Amat, indices, Yvec) 230 376 Adiag = np.sqrt(np.diag(Amat)) 231 377 Anorm = np.outer(Adiag,Adiag) 232 Lam = np.eye(Amat.shape[0])*lam 233 Amatlam = Amat/Anorm 378 Amat = Amat/Anorm 234 379 try: 235 Bmat,Nzero = pinv(Amatlam,xtol) #MoorePenrose inversion (via SVD) & count of zeros 236 if Print: G2fil.G2Print('Found %d SVD zeros'%(Nzero), mode='warn') 380 Bmat,Nzeros = pinv(Amat,xtol) #MoorePenrose inversion (via SVD) & count of zeros 237 381 Bmat = Bmat/Anorm 238 return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[], 239 'SVD0':Nzero,'Converged': ifConverged, 'DelChi2':deltaChi2,'Xvec':Xvec, 'chisq0':chisq00}] 240 except nl.LinAlgError: 241 G2fil.G2Print('ouch #2 linear algebra error in making vcov matrix', mode='error') 242 psing = [] 243 if maxcyc: 244 psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e14)[0]) 245 return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':1,'Xvec':None, 'chisq0':chisq00}] 382 except nl.LinAlgError: # this is unexpected. How did we get this far with a singular matrix? 383 G2fil.G2Print('ouch #5 linear algebra error in making final vcov matrix', mode='error') 384 psing = list(np.where(np.abs(np.diag(nl.qr(Amat)[1])) < 1.e14)[0]) 385 if not len(psing): # make sure at least the worst term is flagged 386 psing = [np.argmin(d)] 387 Amat, indices, Yvec = dropTerms(psing, Amat, indices, Yvec) 388 info = {'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':1,'Xvec':None, 'chisq0':chisqf} 389 info['psing'] = [i for i in range(n) if i not in indices] 390 return [x0,None,info] 391 if Nzeros: # sometimes SVD zeros show up when lam=0 392 d = np.abs(np.diag(nl.qr(Amat)[1])) 393 psing = list(np.where(d < 1.e14)[0]) 394 if len(psing) < Nzeros: 395 psing = list(np.where(d <= sorted(d)[Nzeros1])[0]) 396 if Print: 397 G2fil.G2Print('Found %d SVD zeros w/o Lambda. '+ 398 'Likely problem with variable(s) #%s'%(Nzeros,psing), mode='warn') 399 Amat, indices, Yvec = dropTerms(psing, Amat, indices, Yvec) 400 # expand Bmat by filling with zeros if columns have been dropped 401 psing = [i for i in range(n) if i not in indices] 402 if len(psing): 403 ins = [ji for i,j in enumerate(psing)] 404 Bmat = np.insert(np.insert(Bmat,ins,0,1),ins,0,0) 405 info.update({'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'SVD0':Nzeros, 406 'Converged':ifConverged, 'DelChi2':deltaChi2, 'Xvec':XvecAll, 'chisq0':chisq00}) 407 info['psing'] = [i for i in range(n) if i not in indices] 408 return [x0,Bmat,info] 246 409 247 410 def HessianSVD(func,x0,Hess,args=(),ftol=1.49012e8,xtol=1.e6, maxcyc=0,lamda=3,Print=False,refPlotUpdate=None):
Note: See TracChangeset
for help on using the changeset viewer.