Changeset f168d02 in sasview

Timestamp:

Jul 1, 2008 2:56:32 PM (16 years ago)

Author:

Mathieu Doucet <doucetm@…>

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:

2a92852

Parents:

a17ffdf

Message:

Added slit smearing (still slow)

Location:

pr_inversion

Files:

• c_extensions/Cinvertor.c (modified) (3 diffs)
• c_extensions/invertor.c (modified) (2 diffs)
• c_extensions/invertor.h (modified) (1 diff)
• invertor.py (modified) (5 diffs)
• test/utest_invertor.py (modified) (1 diff)

pr_inversion/c_extensions/Cinvertor.c

@@ -624 +624 @@
 
         iq_value = iq(pars, self->params.d_max, npars, q);
+        return Py_BuildValue("f", iq_value);
+}
+
+const char get_iq_smeared_doc[] =
+        "Function to call to evaluate the scattering intensity.\n"
+        "The scattering intensity is slit-smeared."
+        " @param args: c-parameters, and q\n"
+        " @return: I(q)";
+
+/**
+ * Function to call to evaluate the scattering intensity
+ * The scattering intensity is slit-smeared.
+ * @param args: c-parameters, and q
+ * @return: I(q)
+ */
+static PyObject * get_iq_smeared(Cinvertor *self, PyObject *args) {
+        double *pars;
+        double q, iq_value;
+        PyObject *data_obj;
+        Py_ssize_t npars;
+
+        if (!PyArg_ParseTuple(args, "Od", &data_obj, &q)) return NULL;
+        OUTVECTOR(data_obj,pars,npars);
+
+        iq_value = iq_smeared(pars, self->params.d_max, npars,
+                                                        self->params.slit_height, self->params.slit_width,
+                                                        q, 21);
         return Py_BuildValue("f", iq_value);
 }
@@ -875 +902 @@
                     a[i*nfunc+j] = 1.0/self->params.err[i];
                 } else {
-                        a[i*nfunc+j] = ortho_transformed(self->params.d_max, j+offset, self->params.x[i])/self->params.err[i];
+                        if (self->params.slit_width>0 || self->params.slit_height>0) {
+                                a[i*nfunc+j] = ortho_transformed_smeared(self->params.d_max,
+                                                j+offset, self->params.slit_height, self->params.slit_width,
+                                                self->params.x[i], 21)/self->params.err[i];
+                        } else {
+                                a[i*nfunc+j] = ortho_transformed(self->params.d_max, j+offset, self->params.x[i])/self->params.err[i];
+                        }
                 }
             }
@@ -1009 +1042 @@
                    {"get_nerr", (PyCFunction)get_nerr, METH_VARARGS, get_nerr_doc},
                    {"iq", (PyCFunction)get_iq, METH_VARARGS, get_iq_doc},
+                   {"iq_smeared", (PyCFunction)get_iq_smeared, METH_VARARGS, get_iq_smeared_doc},
                    {"pr", (PyCFunction)get_pr, METH_VARARGS, get_pr_doc},
                    {"get_pr_err", (PyCFunction)get_pr_err, METH_VARARGS, get_pr_err_doc},

pr_inversion/c_extensions/invertor.c

@@ -63 +63 @@
 double ortho_transformed_smeared(double d_max, int n, double height, double width, double q, int npts) {
         double sum, value, y, z;
-        int i, j;
+        int i, j, n_height, n_width;
+        double count_w;
         double fnpts;
         sum = 0.0;
         fnpts = (float)npts-1.0;
 
-        for(i=0; i<npts; i++) {
-                y = -width/2.0+width/fnpts*(float)i;
-                for(j=0; j<npts; j++) {
+        // Check for zero slit size
+        n_height = (height>0) ? npts : 1;
+        n_width  = (width>0)  ? npts : 1;
+
+        count_w = 0.0;
+
+        for(j=0; j<n_height; j++) {
+                if(height>0){
                         z = height/fnpts*(float)j;
+                } else {
+                        z = 0.0;
+                }
+
+                for(i=0; i<n_width; i++) {
+                        if(width>0){
+                                y = -width/2.0+width/fnpts*(float)i;
+                        } else {
+                                y = 0.0;
+                        }
+                        if (((q-y)*(q-y)+z*z)<=0.0) continue;
+                        count_w += 1.0;
                         sum += ortho_transformed(d_max, n, sqrt((q-y)*(q-y)+z*z));
                 }
         }
-
-        return sum/npts/npts/height/width;
+        return sum/count_w;
 }
 
@@ -95 +112 @@
     for (i=0; i<n_c; i++) {
         sum += pars[i] * ortho_transformed(d_max, i+1, q);
+    }
+    return sum;
+}
+
+/**
+ * Scattering intensity calculated from the expansion,
+ * slit-smeared.
+ */
+double iq_smeared(double *pars, double d_max, int n_c, double height, double width, double q, int npts) {
+    double sum = 0.0;
+        int i;
+    for (i=0; i<n_c; i++) {
+        sum += pars[i] * ortho_transformed_smeared(d_max, i+1, height, width, q, npts);
     }
     return sum;

pr_inversion/c_extensions/invertor.h

@@ -55 +55 @@
 double int_pr(double *pars, double d_max, int n_c, int nslice);
 double ortho_transformed_smeared(double d_max, int n, double heigth, double width, double q, int npts);
+double iq_smeared(double *pars, double d_max, int n_c, double height, double width, double q, int npts);
 
 #endif

pr_inversion/invertor.py

@@ -73 +73 @@
     suggested_alpha = 0
     ## Last number of base functions used
-    nfunc = 0
+    nfunc = 10
     ## Last output values
     out = None
@@ -191 +191 @@
         invertor.err = self.err
         invertor.has_bck = self.has_bck
+        invertor.slit_height = self.slit_height
+        invertor.slit_width  = self.slit_width
         
         return invertor
@@ -391 +393 @@
         
         # Construct the a matrix and b vector that represent the problem
+        t_0 = time.time()
         self._get_matrix(nfunc, nq, a, b)
+        #print "elasped: ", time.time()-t_0
              
         # Perform the inversion (least square fit)
@@ -531 +535 @@
         from num_term import Num_terms
         estimator = Num_terms(self.clone())
-        
-        return estimator.num_terms(isquit_func)
+        try:
+            return estimator.num_terms(isquit_func)
+        except:
+            # If we fail, estimate alpha and return the default
+            # number of terms 
+            best_alpha, message, elapsed =self.estimate_alpha(self.nfunc)
+            return self.nfunc, best_alpha, "Could not estimate number of terms"
                     
     def estimate_alpha(self, nfunc):
@@ -572 +581 @@
             # just return the estimate
             if npeaks>1:
-                message = "Your P(r) is not smooth, please check your inversion parameters"
+                #message = "Your P(r) is not smooth, please check your inversion parameters"
+                message = None
                 return pr.suggested_alpha, message, elapsed
             else:

pr_inversion/test/utest_invertor.py

@@ -175 +175 @@
         
         self.assertEqual(self.invertor.test_no_data, None)
+        
+    def test_slitsettings(self):
+        self.invertor.slit_width = 1.0
+        self.assertEqual(self.invertor.slit_width, 1.0)
+        self.invertor.slit_height = 2.0
+        self.assertEqual(self.invertor.slit_height, 2.0)
+        
         
     def test_inversion(self):