Changeset 34f3ad0 in sasview
- Timestamp:
- Apr 27, 2012 9:07:00 AM (13 years ago)
- Branches:
- master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
- Children:
- 8a621ac
- Parents:
- 1f9f3c8a
- Location:
- pr_inversion/src/sans/pr
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
pr_inversion/src/sans/pr/distance_explorer.py
rd9621b17 r34f3ad0 16 16 over that range. 17 17 """ 18 import sys 18 19 19 import numpy20 import math21 import sys22 20 23 21 class Results: … … 40 38 self.bck = [] 41 39 self.d_max = [] 42 ## List of errors found during the last exploration 40 ## List of errors found during the last exploration 43 41 self.errors = [] 42 44 43 45 44 class DistExplorer(object): … … 55 54 56 55 """ 57 self.pr_state 56 self.pr_state = pr_state 58 57 self._default_min = 0.8 * self.pr_state.d_max 59 58 self._default_max = 1.2 * self.pr_state.d_max 60 61 59 62 60 def __call__(self, dmin=None, dmax=None, npts=10): … … 80 78 81 79 # Loop over d_max values 82 for i in range(npts): 83 d = dmin + i * (dmax - dmin) /(npts-1.0)80 for i in range(npts): 81 d = dmin + i * (dmax - dmin) / (npts - 1.0) 84 82 self.pr_state.d_max = d 85 83 try: 86 out, cov = self.pr_state.invert(self.pr_state.nfunc) 84 out, cov = self.pr_state.invert(self.pr_state.nfunc) 87 85 88 86 # Store results … … 100 98 results.pos.append(pos) 101 99 results.pos_err.append(pos_err) 102 results.osc.append(osc) 100 results.osc.append(osc) 103 101 except: 104 102 # This inversion failed, skip this D_max value … … 108 106 109 107 return results 110 -
pr_inversion/src/sans/pr/invertor.py
rea59943 r34f3ad0 20 20 Future work: extend this function to allow topic selection 21 21 """ 22 info_txt = "The inversion approach is based on Moore, J. Appl. Cryst. " 22 info_txt = "The inversion approach is based on Moore, J. Appl. Cryst. " 23 23 info_txt += "(1980) 13, 168-175.\n\n" 24 24 info_txt += "P(r) is set to be equal to an expansion of base functions " … … 62 62 In the following i refers to the ith base function coefficient. 63 63 The matrix has its entries j in its first Npts rows set to 64 A[j][i] = (Fourier transformed base function for point j) 64 A[j][i] = (Fourier transformed base function for point j) 65 65 66 66 We them choose a number of r-points, n_r, to evaluate the second 67 67 derivative of P(r) at. This is used as our regularization term. 68 68 For a vector r of length n_r, the following n_r rows are set to 69 A[j+Npts][i] = (2nd derivative of P(r), d**2(P(r))/d(r)**2, 69 A[j+Npts][i] = (2nd derivative of P(r), d**2(P(r))/d(r)**2, 70 70 evaluated at r[j]) 71 71 … … 101 101 info = {} 102 102 103 104 103 def __init__(self): 105 104 Cinvertor.__init__(self) … … 109 108 restore the state of invertor for pickle 110 109 """ 111 (self.__dict__, self.alpha, self.d_max,112 self.q_min, self.q_max, 110 (self.__dict__, self.alpha, self.d_max, 111 self.q_min, self.q_max, 113 112 self.x, self.y, 114 113 self.err, self.has_bck, … … 120 119 """ 121 120 122 state = (self.__dict__, 121 state = (self.__dict__, 123 122 self.alpha, self.d_max, 124 123 self.q_min, self.q_max, 125 124 self.x, self.y, 126 self.err, self.has_bck, 125 self.err, self.has_bck, 127 126 self.slit_height, self.slit_width, 128 127 ) 129 return (Invertor, tuple(), state, None, None)128 return (Invertor, tuple(), state, None, None) 130 129 131 130 def __setattr__(self, name, value): … … 224 223 225 224 invertor = Invertor() 226 invertor.chi2 = self.chi2 227 invertor.elapsed = self.elapsed 228 invertor.nfunc = self.nfunc 225 invertor.chi2 = self.chi2 226 invertor.elapsed = self.elapsed 227 invertor.nfunc = self.nfunc 229 228 invertor.alpha = self.alpha 230 229 invertor.d_max = self.d_max … … 237 236 invertor.has_bck = self.has_bck 238 237 invertor.slit_height = self.slit_height 239 invertor.slit_width 238 invertor.slit_width = self.slit_width 240 239 241 240 invertor.info = copy.deepcopy(self.info) … … 254 253 In the following i refers to the ith base function coefficient. 255 254 The matrix has its entries j in its first Npts rows set to 256 A[i][j] = (Fourier transformed base function for point j) 255 A[i][j] = (Fourier transformed base function for point j) 257 256 258 257 We them choose a number of r-points, n_r, to evaluate the second … … 271 270 :param nfunc: number of base functions to use. 272 271 :param nr: number of r points to evaluate the 2nd derivative at for the reg. term. 273 :return: c_out, c_cov - the coefficients with covariance matrix 272 :return: c_out, c_cov - the coefficients with covariance matrix 274 273 275 274 """ … … 295 294 The minimization function is set to 296 295 sum_i[ (I_obs(q_i) - I_theo(q_i))/err**2 ] + alpha * reg_term, 297 where the reg_term is given by Svergun: it is the integral of 296 where the reg_term is given by Svergun: it is the integral of 298 297 the square of the first derivative 299 298 of P(r), d(P(r))/dr, integrated over the full range of r. 300 299 301 300 :param nfunc: number of base functions to use. 302 :param nr: number of r points to evaluate the 2nd derivative at 301 :param nr: number of r points to evaluate the 2nd derivative at 303 302 for the reg. term. 304 303 305 :return: c_out, c_cov - the coefficients with covariance matrix 306 307 """ 308 309 #from scipy import optimize 310 #import time 311 304 :return: c_out, c_cov - the coefficients with covariance matrix 305 306 """ 312 307 self.nfunc = nfunc 313 308 # First, check that the current data is valid … … 333 328 334 329 if cov_x is None: 335 cov_x = numpy.ones([nfunc, nfunc])330 cov_x = numpy.ones([nfunc, nfunc]) 336 331 cov_x *= math.fabs(chisqr) 337 332 return out, cov_x … … 343 338 with a set of base functions. 344 339 345 This method is provided as a test. 346 """ 347 #from scipy import optimize 348 340 This method is provided as a test. 341 """ 349 342 # First, check that the current data is valid 350 343 if self.is_valid() <= 0: … … 354 347 p = numpy.ones(nfunc) 355 348 t_0 = time.time() 356 out, cov_x, info, mesg, success= optimize.leastsq(self.pr_residuals, p,349 out, cov_x, _, _, _ = optimize.leastsq(self.pr_residuals, p, 357 350 full_output=1) 358 351 … … 371 364 372 365 def pr_err(self, c, c_cov, r): 373 """ 366 """ 374 367 Returns the value of P(r) for a given r, and base function 375 368 coefficients, with error. … … 403 396 In the following i refers to the ith base function coefficient. 404 397 The matrix has its entries j in its first Npts rows set to 405 A[i][j] = (Fourier transformed base function for point j) 398 A[i][j] = (Fourier transformed base function for point j) 406 399 407 400 We them choose a number of r-points, n_r, to evaluate the second … … 449 442 nfunc += 1 450 443 451 a = numpy.zeros([npts +nq, nfunc])452 b = numpy.zeros(npts +nq)444 a = numpy.zeros([npts + nq, nfunc]) 445 b = numpy.zeros(npts + nq) 453 446 err = numpy.zeros([nfunc, nfunc]) 454 447 … … 461 454 462 455 # Perform the inversion (least square fit) 463 c, chi2, rank, n= lstsq(a, b)456 c, chi2, _, _ = lstsq(a, b) 464 457 # Sanity check 465 458 try: … … 469 462 self.chi2 = chi2 470 463 471 inv_cov = numpy.zeros([nfunc, nfunc])464 inv_cov = numpy.zeros([nfunc, nfunc]) 472 465 # Get the covariance matrix, defined as inv_cov = a_transposed * a 473 466 self._get_invcov_matrix(nfunc, nr, a, inv_cov) … … 477 470 478 471 if math.fabs(self.alpha) > 0: 479 new_alpha = sum_sig /(sum_reg/self.alpha)472 new_alpha = sum_sig / (sum_reg / self.alpha) 480 473 else: 481 474 new_alpha = 0.0 … … 484 477 try: 485 478 cov = numpy.linalg.pinv(inv_cov) 486 err = math.fabs(chi2 /float(npts-nfunc)) * cov479 err = math.fabs(chi2 / float(npts - nfunc)) * cov 487 480 except: 488 481 # We were not able to estimate the errors … … 509 502 self.cov = err_0 510 503 504 # Store computation time 505 self.elapsed = time.time() - t_0 506 511 507 return self.out, self.cov 512 508 … … 516 512 number of terms 517 513 518 :param isquit_func: reference to thread function to call to 514 :param isquit_func: reference to thread function to call to 519 515 check whether the computation needs to 520 516 be stopped. … … 529 525 except: 530 526 # If we fail, estimate alpha and return the default 531 # number of terms 532 best_alpha, message, elapsed =self.estimate_alpha(self.nfunc)527 # number of terms 528 best_alpha, _, _ = self.estimate_alpha(self.nfunc) 533 529 return self.nfunc, best_alpha, "Could not estimate number of terms" 534 530 … … 547 543 """ 548 544 #import time 549 try: 545 try: 550 546 pr = self.clone() 551 547 … … 560 556 561 557 # Perform inversion to find the largest alpha 562 out, cov= pr.invert(nfunc)558 out, _ = pr.invert(nfunc) 563 559 elapsed = time.time() - starttime 564 560 initial_alpha = pr.alpha … … 567 563 # Try the inversion with the estimated alpha 568 564 pr.alpha = pr.suggested_alpha 569 out, cov= pr.invert(nfunc)565 out, _ = pr.invert(nfunc) 570 566 571 567 npeaks = pr.get_peaks(out) … … 573 569 # just return the estimate 574 570 if npeaks > 1: 575 #message = "Your P(r) is not smooth, 571 #message = "Your P(r) is not smooth, 576 572 #please check your inversion parameters" 577 573 message = None … … 587 583 for i in range(10): 588 584 pr.alpha = (0.33)**(i+1) * alpha 589 out, cov= pr.invert(nfunc)585 out, _ = pr.invert(nfunc) 590 586 591 587 peaks = pr.get_peaks(out) … … 599 595 # just return that 600 596 if not found and initial_peaks == 1 and \ 601 initial_alpha <best_alpha:597 initial_alpha < best_alpha: 602 598 best_alpha = initial_alpha 603 599 … … 608 604 message = None 609 605 elif best_alpha >= 0.5 * pr.suggested_alpha: 610 # best alpha is too big, return a 606 # best alpha is too big, return a 611 607 # reasonable value 612 608 message = "The estimated alpha for your system is too " 613 mess sage += "large. "609 message += "large. " 614 610 message += "Try increasing your maximum distance." 615 611 … … 620 616 return 0, message, elapsed 621 617 622 623 618 def to_file(self, path, npts=100): 624 619 """ … … 660 655 file.close() 661 656 662 663 657 def from_file(self, path): 664 658 """ … … 676 670 fd = open(path, 'r') 677 671 678 buff 679 lines 672 buff = fd.read() 673 lines = buff.split('\n') 680 674 for line in lines: 681 675 if line.startswith('#d_max='): … … 736 730 self.out[i] = float(toks2[0]) 737 731 738 self.cov[i][i] = float(toks2[1]) 732 self.cov[i][i] = float(toks2[1]) 739 733 740 734 except: … … 742 736 raise RuntimeError, msg 743 737 else: 744 msg = "Invertor.from_file: '%s' is not a file" % str(path) 738 msg = "Invertor.from_file: '%s' is not a file" % str(path) 745 739 raise RuntimeError, msg 746 747 748 if __name__ == "__main__":749 o = Invertor()750 751 752 753 754 -
pr_inversion/src/sans/pr/num_term.py
r3fd8390 r34f3ad0 1 2 #import unittest3 1 import math 4 2 import numpy 5 #import sys6 #import string7 3 import copy 8 4 from sans.pr.invertor import Invertor 9 5 10 class Num_terms(): 6 7 class Num_terms(): 11 8 """ 12 9 """ … … 62 59 def get0_out(self): 63 60 """ 64 """ 61 """ 65 62 inver = self.invertor 66 63 self.osc_list = [] … … 70 67 if self.isquit_func != None: 71 68 self.isquit_func() 72 best_alpha, message, elapsed= inver.estimate_alpha(k)69 best_alpha, message, _ = inver.estimate_alpha(k) 73 70 inver.alpha = best_alpha 74 71 inver.out, inver.cov = inver.lstsq(k) … … 82 79 self.mess_list.append(message) 83 80 84 #print "osc", self.osc_list85 #print "err", self.err_list86 #print "alpha", self.alpha_list87 #new_ls = []88 81 new_osc1 = [] 89 82 new_osc2 = [] … … 101 94 if self.err_list[i] < 0.8 and self.err_list[i] >= 0.7: 102 95 new_osc3.append(self.osc_list[i]) 103 f alg7 = True96 flag7 = True 104 97 105 98 if flag9 == True: … … 110 103 self.dataset = new_osc3 111 104 112 #print "dataset", self.dataset113 105 return self.dataset 114 106 … … 148 140 try: 149 141 self.isquit_func = isquit_func 150 #self.nterm_max = len(self.invertor.x)151 #self.nterm_max = 32152 142 nts = self.compare_err() 153 #print "nts", nts154 143 div = len(nts) 155 144 tem = float(div)/2.0 … … 163 152 return self.nterm_min, self.alpha_list[10], self.mess_list[10] 164 153 154 165 155 #For testing 166 156 def load(path): 167 #import numpy, math, sys168 157 # Read the data from the data file 169 158 data_x = numpy.zeros(0) … … 187 176 scale = 0.05 * math.sqrt(test_y) 188 177 #scale = 0.05/math.sqrt(y) 189 min_err = 0.01 *y178 min_err = 0.01 * y 190 179 err = scale * math.sqrt(test_y) + min_err 191 180 #err = 0 … … 204 193 x, y, erro = load("test/Cyl_A_D102.txt") 205 194 i.d_max = 102.0 206 i.nfunc = 10 195 i.nfunc = 10 207 196 #i.q_max = 0.4 208 197 #i.q_min = 0.07 209 i.x 210 i.y 211 i.err = erro 198 i.x = x 199 i.y = y 200 i.err = erro 212 201 #i.out, i.cov = i.lstsq(10) 213 202 # Testing estimator 214 203 est = Num_terms(i) 215 204 print est.num_terms() 216 217 218 219 220 221
Note: See TracChangeset
for help on using the changeset viewer.