-                      r79492222
+                      r3350ad6
+# pylint: disable=invalid-name
 """
 Module to perform P(r) inversion.
 …
 import os
 import re
+import logging
 from numpy.linalg import lstsq
 from scipy import optimize
 …
     Future work: extend this function to allow topic selection
     """
     info_txt  = "The inversion approach is based on Moore, J. Appl. Cryst. "
+    info_txt = "The inversion approach is based on Moore, J. Appl. Cryst. "
     info_txt += "(1980) 13, 168-175.\n\n"
     info_txt += "P(r) is set to be equal to an expansion of base functions "
 …
     info_txt += "   - Maximum distance: the maximum distance between any "
     info_txt += "two points in the system.\n"
     return info_txt
 class Invertor(Cinvertor):
     """
     Invertor class to perform P(r) inversion
     The problem is solved by posing the problem as  Ax = b,
     where x is the set of coefficients we are looking for.
     Npts is the number of points.
     In the following i refers to the ith base function coefficient.
     The matrix has its entries j in its first Npts rows set to ::
         A[j][i] = (Fourier transformed base function for point j)
     We them choose a number of r-points, n_r, to evaluate the second
     derivative of P(r) at. This is used as our regularization term.
 …
         A[j+Npts][i] = (2nd derivative of P(r), d**2(P(r))/d(r)**2,
         evaluated at r[j])
     The vector b has its first Npts entries set to ::
         b[j] = (I(q) observed for point j)
     The following n_r entries are set to zero.
     The result is found by using scipy.linalg.basic.lstsq to invert
     the matrix and find the coefficients x.
     Methods inherited from Cinvertor:
 …
     """
     ## Chisqr of the last computation
     chi2  = 0
+    chi2 = 0
     ## Time elapsed for last computation
     elapsed = 0
     ## Alpha to get the reg term the same size as the signal
     suggested_alpha = 0
+    alpha = 0
     ## Last number of base functions used
     nfunc = 10
 …
     ## Information dictionary for application use
     info = {}
     def __init__(self):
         Cinvertor.__init__(self)
+        self.slit_height = None
+        self.slit_width = None
+        self.err = None
+        self.d_max = None
+        self.q_min = self.set_qmin(-1.0)
+        self.q_max = self.set_qmax(-1.0)
+        self.x = None
+        self.y = None
+        self.has_bck = False
     def __setstate__(self, state):
         """
 …
          self.err, self.has_bck,
          self.slit_height, self.slit_width) = state
     def __reduce_ex__(self, proto):
         """
 …
                  self.err, self.has_bck,
                  self.slit_height, self.slit_width,
+                 )
+                )
         return (Invertor, tuple(), state, None, None)
     def __setattr__(self, name, value):
         """
 …
             else:
                 raise ValueError, "Invertor: has_bck can only be True or False"
         return Cinvertor.__setattr__(self, name, value)
     def __getattr__(self, name):
         """
 …
             return self.__dict__[name]
         return None
     def clone(self):
         """
 …
         """
         #import copy
         invertor = Invertor()
         invertor.chi2    = self.chi2
+        invertor.chi2 = self.chi2
         invertor.elapsed = self.elapsed
         invertor.nfunc   = self.nfunc
         invertor.alpha   = self.alpha
         invertor.d_max   = self.d_max
         invertor.q_min   = self.q_min
         invertor.q_max   = self.q_max
+        invertor.nfunc = self.nfunc
+        invertor.alpha = self.alpha
+        invertor.d_max = self.d_max
+        invertor.q_min = self.q_min
+        invertor.q_max = self.q_max
         invertor.x = self.x
         invertor.y = self.y
 …
         invertor.slit_height = self.slit_height
         invertor.slit_width = self.slit_width
         invertor.info = copy.deepcopy(self.info)
         return invertor
     def invert(self, nfunc=10, nr=20):
         """
         Perform inversion to P(r)
         The problem is solved by posing the problem as  Ax = b,
         where x is the set of coefficients we are looking for.
         Npts is the number of points.
         In the following i refers to the ith base function coefficient.
         The matrix has its entries j in its first Npts rows set to ::
             A[i][j] = (Fourier transformed base function for point j)
         We them choose a number of r-points, n_r, to evaluate the second
         derivative of P(r) at. This is used as our regularization term.
 …
             A[i+Npts][j] = (2nd derivative of P(r), d**2(P(r))/d(r)**2, evaluated at r[j])
         The vector b has its first Npts entries set to ::
             b[j] = (I(q) observed for point j)
         The following n_r entries are set to zero.
         The result is found by using scipy.linalg.basic.lstsq to invert
         the matrix and find the coefficients x.
         :param nfunc: number of base functions to use.
         :param nr: number of r points to evaluate the 2nd derivative at for the reg. term.
 …
         self.background = 0.0
         return self.lstsq(nfunc, nr=nr)
     def iq(self, out, q):
         """
         Function to call to evaluate the scattering intensity
         :param args: c-parameters, and q
         :return: I(q)
         """
         return Cinvertor.iq(self, out, q) + self.background
     def invert_optimize(self, nfunc=10, nr=20):
         """
         Slower version of the P(r) inversion that uses scipy.optimize.leastsq.
         This probably produce more reliable results, but is much slower.
         The minimization function is set to
 …
         the square of the first derivative
         of P(r), d(P(r))/dr, integrated over the full range of r.
         :param nfunc: number of base functions to use.
         :param nr: number of r points to evaluate the 2nd derivative at
             for the reg. term.
         :return: c_out, c_cov - the coefficients with covariance matrix
         """
         self.nfunc = nfunc
 …
             msg = "Invertor.invert: Data array are of different length"
             raise RuntimeError, msg
         p = numpy.ones(nfunc)
         t_0 = time.time()
+        out, cov_x, _, _, _ = optimize.leastsq(self.residuals,
+                                                            p, full_output=1)
+        out, cov_x, _, _, _ = optimize.leastsq(self.residuals, p, full_output=1)
         # Compute chi^2
         res = self.residuals(out)
 …
         for i in range(len(res)):
             chisqr += res[i]
         self.chi2 = chisqr
         # Store computation time
         self.elapsed = time.time() - t_0
         if cov_x is None:
             cov_x = numpy.ones([nfunc, nfunc])
             cov_x *= math.fabs(chisqr)
         return out, cov_x
     def pr_fit(self, nfunc=5):
         """
 …
         is set to some P(r) distribution that we are trying to reproduce
         with a set of base functions.
         This method is provided as a test.
         """
 …
             msg = "Invertor.invert: Data arrays are of different length"
             raise RuntimeError, msg
         p = numpy.ones(nfunc)
         t_0 = time.time()
+        out, cov_x, _, _, _ = optimize.leastsq(self.pr_residuals, p,
+                                                            full_output=1)
+        out, cov_x, _, _, _ = optimize.leastsq(self.pr_residuals, p, full_output=1)
         # Compute chi^2
         res = self.pr_residuals(out)
 …
         for i in range(len(res)):
             chisqr += res[i]
         self.chisqr = chisqr
         # Store computation time
         self.elapsed = time.time() - t_0
         return out, cov_x
     def pr_err(self, c, c_cov, r):
         """
         Returns the value of P(r) for a given r, and base function
         coefficients, with error.
         :param c: base function coefficients
         :param c_cov: covariance matrice of the base function coefficients
         :param r: r-value to evaluate P(r) at
         :return: P(r)
         """
         return self.get_pr_err(c, c_cov, r)
     def _accept_q(self, q):
         """
 …
             return False
         return True
     def lstsq(self, nfunc=5, nr=20):
         """
         The problem is solved by posing the problem as  Ax = b,
         where x is the set of coefficients we are looking for.
         Npts is the number of points.
         In the following i refers to the ith base function coefficient.
         The matrix has its entries j in its first Npts rows set to ::
             A[i][j] = (Fourier transformed base function for point j)
         We them choose a number of r-points, n_r, to evaluate the second
         derivative of P(r) at. This is used as our regularization term.
 …
             A[i+Npts][j] = (2nd derivative of P(r), d**2(P(r))/d(r)**2,
             evaluated at r[j])
         The vector b has its first Npts entries set to ::
             b[j] = (I(q) observed for point j)
         The following n_r entries are set to zero.
         The result is found by using scipy.linalg.basic.lstsq to invert
         the matrix and find the coefficients x.
         :param nfunc: number of base functions to use.
         :param nr: number of r points to evaluate the 2nd derivative at for the reg. term.
 …
         # blah = numpy.ascontiguousarray(blah_original)
         # ... before passing it to C
         if self.is_valid() < 0:
             msg = "Invertor: invalid data; incompatible data lengths."
             raise RuntimeError, msg
         self.nfunc = nfunc
         # a -- An M x N matrix.
         # b -- An M x nrhs matrix or M vector.
         npts = len(self.x)
         nq   = nr
+        nq = nr
         sqrt_alpha = math.sqrt(math.fabs(self.alpha))
         if sqrt_alpha < 0.0:
 …
         b = numpy.zeros(npts + nq)
         err = numpy.zeros([nfunc, nfunc])
         # Construct the a matrix and b vector that represent the problem
         t_0 = time.time()
 …
         except:
             raise RuntimeError, "Invertor: could not invert I(Q)\n  %s" % sys.exc_value
         # Perform the inversion (least square fit)
         c, chi2, _, _ = lstsq(a, b)
 …
             chi2 = -1.0
         self.chi2 = chi2
         inv_cov = numpy.zeros([nfunc, nfunc])
         # Get the covariance matrix, defined as inv_cov = a_transposed * a
         self._get_invcov_matrix(nfunc, nr, a, inv_cov)
         # Compute the reg term size for the output
         sum_sig, sum_reg = self._get_reg_size(nfunc, nr, a)
         if math.fabs(self.alpha) > 0:
             new_alpha = sum_sig / (sum_reg / self.alpha)
 …
             new_alpha = 0.0
         self.suggested_alpha = new_alpha
         try:
             cov = numpy.linalg.pinv(inv_cov)
 …
             # We were not able to estimate the errors
             # Return an empty error matrix
             pass
+            logging.error(sys.exc_value)
         # Keep a copy of the last output
         if self.has_bck == False:
 …
         else:
             self.background = c[0]
             err_0 = numpy.zeros([nfunc, nfunc])
             c_0 = numpy.zeros(nfunc)
             for i in range(nfunc_0):
                 c_0[i] = c[i+1]
+                c_0[i] = c[i + 1]
                 for j in range(nfunc_0):
                     err_0[i][j] = err[i+1][j+1]
+                    err_0[i][j] = err[i + 1][j + 1]
             self.out = c_0
             self.cov = err_0
         # Store computation time
         self.elapsed = time.time() - t_0
         return self.out, self.cov
     def estimate_numterms(self, isquit_func=None):
         """
         Returns a reasonable guess for the
         number of terms
         :param isquit_func:
           reference to thread function to call to check whether the computation needs to
           be stopped.
         :return: number of terms, alpha, message
         """
         from num_term import Num_terms
 …
             best_alpha, _, _ = self.estimate_alpha(self.nfunc)
             return self.nfunc, best_alpha, "Could not estimate number of terms"
     def estimate_alpha(self, nfunc):
         """
         Returns a reasonable guess for the
         regularization constant alpha
         :param nfunc: number of terms to use in the expansion.
         :return: alpha, message, elapsed
         where alpha is the estimate for alpha,
         message is a message for the user,
 …
         try:
             pr = self.clone()
             # T_0 for computation time
             starttime = time.time()
             elapsed = 0
             # If the current alpha is zero, try
             # another value
             if pr.alpha <= 0:
                 pr.alpha = 0.0001
             # Perform inversion to find the largest alpha
             out, _ = pr.invert(nfunc)
 …
             initial_alpha = pr.alpha
             initial_peaks = pr.get_peaks(out)
             # Try the inversion with the estimated alpha
             pr.alpha = pr.suggested_alpha
             out, _ = pr.invert(nfunc)
             npeaks = pr.get_peaks(out)
             # if more than one peak to start with
 …
                 return pr.suggested_alpha, message, elapsed
             else:
                 # Look at smaller values
                 # We assume that for the suggested alpha, we have 1 peak
 …
                 found = False
                 for i in range(10):
                     pr.alpha = (0.33)**(i+1) * alpha
+                    pr.alpha = (0.33) ** (i + 1) * alpha
                     out, _ = pr.invert(nfunc)
                     peaks = pr.get_peaks(out)
                     if peaks > 1:
 …
                         break
                     best_alpha = pr.alpha
                 # If we didn't find a turning point for alpha and
                 # the initial alpha already had only one peak,
 …
                     initial_alpha < best_alpha:
                     best_alpha = initial_alpha
                 # Check whether the size makes sense
                 message = ''
                 if not found:
                     message = None
 …
                     # best alpha is too big, return a
                     # reasonable value
                     message  = "The estimated alpha for your system is too "
+                    message = "The estimated alpha for your system is too "
                     message += "large. "
                     message += "Try increasing your maximum distance."
                 return best_alpha, message, elapsed
         except:
             message = "Invertor.estimate_alpha: %s" % sys.exc_value
             return 0, message, elapsed
     def to_file(self, path, npts=100):
         """
         Save the state to a file that will be readable
         by SliceView.
         :param path: path of the file to write
         :param npts: number of P(r) points to be written
         """
         file = open(path, 'w')
 …
                 for i in range(len(self.out)):
                     file.write("#C_%i=%s+-%s\n" % (i, str(self.out[i]),
                                                     str(self.cov[i][i])))
+                                                   str(self.cov[i][i])))
         file.write("<r>  <Pr>  <dPr>\n")
         r = numpy.arange(0.0, self.d_max, self.d_max/npts)
+        r = numpy.arange(0.0, self.d_max, self.d_max / npts)
         for r_i in r:
             (value, err) = self.pr_err(self.out, self.cov, r_i)
             file.write("%g  %g  %g\n" % (r_i, value, err))
         file.close()
     def from_file(self, path):
         """
 …
         to be able to generate P(r) from a set of
         parameters.
         :param path: path of the file to load
         """
         #import os
 …
             try:
                 fd = open(path, 'r')
                 buff = fd.read()
                 lines = buff.split('\n')
 …
                         else:
                             self.has_bck = False
                     # Now read in the parameters
                     elif line.startswith('#C_'):
 …
                         i = int(m.group(1))
                         self.out[i] = float(toks2[0])
                         self.cov[i][i] = float(toks2[1])
             except:
                 msg = "Invertor.from_file: corrupted file\n%s" % sys.exc_value

Note: See TracChangeset for help on using the changeset viewer.

SasView

Changeset 3350ad6 in sasview

Legend:

src/sas/pr/invertor.py

Download in other formats: