Changeset b8080e1 in sasview for src/sas/sascalc/calculator

Timestamp:

Aug 29, 2018 10:01:23 AM (6 years ago)

Author:

Piotr Rozyczko <rozyczko@…>

Branches:

ESS_GUI, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc

Children:

9463ca2

Parents:

ce30949

git-author:

Piotr Rozyczko <rozyczko@…> (08/29/18 09:59:56)

git-committer:

Piotr Rozyczko <rozyczko@…> (08/29/18 10:01:23)

Message:

cherry picking sascalc changes from master SASVIEW-996
minor unit test fixes

Location:

src/sas/sascalc/calculator

Files:

• c_extensions/librefl.c (modified) (5 diffs)
• c_extensions/sld2i.c (modified) (8 diffs)
• c_extensions/sld2i.h (modified) (2 diffs)
• c_extensions/sld2i_module.c (modified) (5 diffs)
• instrument.py (modified) (1 diff)
• resolution_calculator.py (modified) (3 diffs)
• sas_gen.py (modified) (9 diffs)

src/sas/sascalc/calculator/c_extensions/librefl.c

@@ -7 +7 @@
 #include <stdio.h>
 #include <stdlib.h>
-#if defined(_MSC_VER)
+#if defined _MSC_VER  || defined __TINYCC__
 #define NEED_ERF
 #endif
@@ -21 +21 @@
 
 
-#ifdef _WIN32
+#ifdef __TINYCC__
+# ifdef isnan
+#   undef isnan
+# endif
+# ifdef isfinite
+#   undef isfinite
+# endif
+# define isnan(x) (x != x)
+# define isfinite(x) (x != INFINITY && x != -INFINITY)
+#elif defined _WIN32
 # include <float.h>
 # if !defined __MINGW32__ || defined __NO_ISOCEXT
@@ -30 +39 @@
 #   define isinf(x) (!_finite(x) && !_isnan(x))
 #  endif
-#  ifndef finite
-#   define finite(x) _finite(x)
+#  ifndef isfinite
+#   define isfinite(x) _finite(x)
 #  endif
 # endif
@@ -84 +93 @@
 double erf(double x)
 {
-    if (!finite(x)) {
+    if (!isfinite(x)) {
         if (isnan(x)) return x;      /* erf(NaN)   = NaN   */
         return (x>0 ? 1.0 : -1.0);   /* erf(+-inf) = +-1.0 */
@@ -94 +103 @@
 double erfc(double x)
 {
-    if (!finite(x)) {
+    if (!isfinite(x)) {
         if (isnan(x)) return x;      /* erfc(NaN)   = NaN      */
         return (x>0 ? 0.0 : 2.0);    /* erfc(+-inf) = 0.0, 2.0 */

src/sas/sascalc/calculator/c_extensions/sld2i.c

@@ -25 +25 @@
  * @param s_theta: angle (from x-axis) of the up spin in degree
  */
-void initGenI(GenI* this, int npix, double* x, double* y, double* z, double* sldn,
+void initGenI(GenI* this, int is_avg, int npix, double* x, double* y, double* z, double* sldn,
                         double* mx, double* my, double* mz, double* voli,
                         double in_spin, double out_spin,
                         double s_theta) {
+        this->is_avg = is_avg;
         this->n_pix = npix;
         this->x_val = x;
@@ -63 +64 @@
         Cplx temp_fi;
 
+        double count = 0.0;
         int i, j;
-
-        double count = 0.0;
-        //check if this computation is for averaging
 
         cassign(&iqr, 0.0, 0.0);
@@ -78 +77 @@
 
         // Loop over q-values and multiply apply matrix
+
+        //printf("npoints: %d, npix: %d\n", npoints, this->n_pix);
         for(i=0; i<npoints; i++){
                 //I_out[i] = 0.0;
@@ -86 +87 @@
                 //printf("i: %d\n", i);
                 //q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]); // + qz[i]*qz[i]);
+
                 for(j=0; j<this->n_pix; j++){
                         if (this->sldn_val[j]!=0.0
@@ -143 +145 @@
                 I_out[i] *= (1.0E+8 / count); //in cm (unit) / number; //to be multiplied by vol_pix
         }
-        //printf ("count = %d %g %g %g %g\n", count, sldn_val[0],mx_val[0], my_val[0], mz_val[0]);
+        //printf("count = %d %g %g %g %g\n", count, this->sldn_val[0],this->mx_val[0], this->my_val[0], this->mz_val[0]);
 }
 /**
@@ -154 +156 @@
         // Assumes that q doesn't have qz component and sld_n is all real
         //double Pi = 4.0*atan(1.0);
-        int is_sym = this->n_pix < 0;
         double qr = 0.0;
         double sumj;
         double sld_j = 0.0;
         double count = 0.0;
-        int n_pix = is_sym ? -this->n_pix : this->n_pix;
+        int i, j, k;
+
         //Assume that pixel volumes are given in vol_pix in A^3 unit
         // Loop over q-values and multiply apply matrix
-    int i, j, k;
         for(i=0; i<npoints; i++){
                 sumj =0.0;
-                for(j=0; j<n_pix; j++){
+                for(j=0; j<this->n_pix; j++){
                         //Isotropic: Assumes all slds are real (no magnetic)
                         //Also assumes there is no polarization: No dependency on spin
-                        if (is_sym == 1){
+                        if (this->is_avg == 1){
                                 // approximation for a spherical symmetric particle
                                 qr = sqrt(this->x_val[j]*this->x_val[j]+this->y_val[j]*this->y_val[j]+this->z_val[j]*this->z_val[j])*q[i];
@@ -182 +183 @@
                                 //full calculation
                                 //pragma omp parallel for
-                                for(k=0; k<n_pix; k++){
+                                for(k=0; k<this->n_pix; k++){
                                         sld_j =  this->sldn_val[j] * this->sldn_val[k] * this->vol_pix[j] * this->vol_pix[k];
                                         qr = (this->x_val[j]-this->x_val[k])*(this->x_val[j]-this->x_val[k])+
@@ -201 +202 @@
                 }
                 I_out[i] = sumj;
-                if (is_sym == 1){
+                if (this->is_avg == 1) {
                         I_out[i] *= sumj;
                 }
                 I_out[i] *= (1.0E+8 / count); //in cm (unit) / number; //to be multiplied by vol_pix
         }
-        //printf ("count = %d %g %g %g %g\n", count, sldn_val[0],mx_val[0], my_val[0], mz_val[0]);
+        //printf("count = %d %g %g %g %g\n", count, sldn_val[0],mx_val[0], my_val[0], mz_val[0]);
 }

src/sas/sascalc/calculator/c_extensions/sld2i.h

@@ -10 +10 @@
 typedef struct {
         // vectors
+        int is_avg;
         int n_pix;
         double* x_val;
@@ -26 +27 @@
 
 // Constructor
-void initGenI(GenI*, int npix, double* x, double* y, double* z, double* sldn,
-                double* mx, double* my, double* mz, double* voli,
+void initGenI(GenI*, int is_avg, int npix, double* x, double* y, double* z,
+                double* sldn, double* mx, double* my, double* mz, double* voli,
                 double in_spin, double out_spin,
                 double s_theta);

src/sas/sascalc/calculator/c_extensions/sld2i_module.c

@@ -35 +35 @@
 void
 del_sld2i(PyObject *obj){
+#if PY_MAJOR_VERSION < 3
+        GenI* sld2i = (GenI *)obj;
+#else
         GenI* sld2i = (GenI *)(PyCapsule_GetPointer(obj, "GenI"));
+#endif
         PyMem_Free((void *)sld2i);
 }
@@ -51 +55 @@
         PyObject *mz_val_obj;
         PyObject *vol_pix_obj;
-        Py_ssize_t n_x;
-        //PyObject rlimit_obj;
-        //PyObject npoints_obj;
-        //PyObject nrbins_obj;
-        //PyObject nphibins_obj;
-        int n_pix;
+        Py_ssize_t n_x, n_y, n_z, n_sld, n_mx, n_my, n_mz, n_vol_pix;
+        int is_avg;
         double* x_val;
         double* y_val;
@@ -68 +68 @@
         double outspin;
         double stheta;
-        GenI *sld2i;
-
-        if (!PyArg_ParseTuple(args, "iOOOOOOOOddd", &n_pix, &x_val_obj, &y_val_obj, &z_val_obj, &sldn_val_obj, &mx_val_obj, &my_val_obj, &mz_val_obj, &vol_pix_obj, &inspin, &outspin, &stheta)) return NULL;
-        OUTVECTOR(x_val_obj, x_val, n_x);
-        OUTVECTOR(y_val_obj, y_val, n_x);
-        OUTVECTOR(z_val_obj, z_val, n_x);
-        OUTVECTOR(sldn_val_obj, sldn_val, n_x);
-        OUTVECTOR(mx_val_obj, mx_val, n_x);
-        OUTVECTOR(my_val_obj, my_val, n_x);
-        OUTVECTOR(mz_val_obj, mz_val, n_x);
-        OUTVECTOR(vol_pix_obj, vol_pix, n_x);
-        sld2i =  PyMem_Malloc(sizeof(GenI));
+        PyObject *obj;
+        GenI* sld2i;
+
+        //printf("new GenI\n");
+        if (!PyArg_ParseTuple(args, "iOOOOOOOOddd", &is_avg, &x_val_obj, &y_val_obj, &z_val_obj, &sldn_val_obj, &mx_val_obj, &my_val_obj, &mz_val_obj, &vol_pix_obj, &inspin, &outspin, &stheta)) return NULL;
+        INVECTOR(x_val_obj, x_val, n_x);
+        INVECTOR(y_val_obj, y_val, n_y);
+        INVECTOR(z_val_obj, z_val, n_z);
+        INVECTOR(sldn_val_obj, sldn_val, n_sld);
+        INVECTOR(mx_val_obj, mx_val, n_mx);
+        INVECTOR(my_val_obj, my_val, n_my);
+        INVECTOR(mz_val_obj, mz_val, n_mz);
+        INVECTOR(vol_pix_obj, vol_pix, n_vol_pix);
+        sld2i = PyMem_Malloc(sizeof(GenI));
+        //printf("sldi:%p\n", sld2i);
         if (sld2i != NULL) {
-                initGenI(sld2i, n_pix,x_val,y_val,z_val,sldn_val,mx_val,my_val,mz_val,vol_pix,inspin,outspin,stheta);
+                initGenI(sld2i,is_avg,(int)n_x,x_val,y_val,z_val,sldn_val,mx_val,my_val,mz_val,vol_pix,inspin,outspin,stheta);
         }
-        return PyCapsule_New(sld2i, "GenI", del_sld2i);
+        obj = PyCapsule_New(sld2i, "GenI", del_sld2i);
+        //printf("constructed %p\n", obj);
+        return obj;
 }
 
@@ -90 +95 @@
  */
 PyObject * genicom_inputXY(PyObject *self, PyObject *args) {
-        int npoints;
+        PyObject *gen_obj;
         PyObject *qx_obj;
+        PyObject *qy_obj;
+        PyObject *I_out_obj;
+        Py_ssize_t n_qx, n_qy, n_out;
         double *qx;
-        PyObject *qy_obj;
         double *qy;
-        PyObject *I_out_obj;
-        Py_ssize_t n_out;
         double *I_out;
-        PyObject *gen_obj;
-        GenI *sld2i;
-
-        if (!PyArg_ParseTuple(args, "OiOOO",  &gen_obj, &npoints, &qx_obj, &qy_obj, &I_out_obj)) return NULL;
-        OUTVECTOR(qx_obj, qx, n_out);
-        OUTVECTOR(qy_obj, qy, n_out);
+        GenI* sld2i;
+
+        //printf("in genicom_inputXY\n");
+        if (!PyArg_ParseTuple(args, "OOOO",  &gen_obj, &qx_obj, &qy_obj, &I_out_obj)) return NULL;
+        sld2i = (GenI *)PyCapsule_GetPointer(gen_obj, "GenI");
+        INVECTOR(qx_obj, qx, n_qx);
+        INVECTOR(qy_obj, qy, n_qy);
         OUTVECTOR(I_out_obj, I_out, n_out);
+        //printf("qx, qy, I_out: %d %d %d, %d %d %d\n", qx, qy, I_out, n_qx, n_qy, n_out);
 
         // Sanity check
-        //if(n_in!=n_out) return Py_BuildValue("i",-1);
-
-        // Set the array pointers
-        sld2i = (GenI *)PyCapsule_GetPointer(gen_obj, "GenI");
-
-        genicomXY(sld2i, npoints, qx, qy, I_out);
+        //if(n_q!=n_out) return Py_BuildValue("i",-1);
+
+        genicomXY(sld2i, (int)n_qx, qx, qy, I_out);
+        //printf("done calc\n");
         //return PyCObject_FromVoidPtr(s, del_genicom);
         return Py_BuildValue("i",1);
@@ -121 +126 @@
  */
 PyObject * genicom_input(PyObject *self, PyObject *args) {
-        int npoints;
+        PyObject *gen_obj;
         PyObject *q_obj;
+        PyObject *I_out_obj;
+        Py_ssize_t n_q, n_out;
         double *q;
-        PyObject *I_out_obj;
-        Py_ssize_t n_out;
         double *I_out;
-        PyObject *gen_obj;
         GenI *sld2i;
 
-        if (!PyArg_ParseTuple(args, "OiOO",  &gen_obj, &npoints, &q_obj, &I_out_obj)) return NULL;
-        OUTVECTOR(q_obj, q, n_out);
+        if (!PyArg_ParseTuple(args, "OOO",  &gen_obj, &q_obj, &I_out_obj)) return NULL;
+        sld2i = (GenI *)PyCapsule_GetPointer(gen_obj, "GenI");
+        INVECTOR(q_obj, q, n_q);
         OUTVECTOR(I_out_obj, I_out, n_out);
 
         // Sanity check
-        //if(n_in!=n_out) return Py_BuildValue("i",-1);
-
-        // Set the array pointers
-        sld2i = (GenI *)PyCapsule_GetPointer(gen_obj, "GenI");
-
-        genicom(sld2i, npoints, q, I_out);
-        //return PyCObject_FromVoidPtr(s, del_genicom);
+        //if (n_q!=n_out) return Py_BuildValue("i",-1);
+
+        genicom(sld2i, (int)n_q, q, I_out);
         return Py_BuildValue("i",1);
 }

src/sas/sascalc/calculator/instrument.py

@@ -314 +314 @@
         """
         To plot the wavelength spactrum
-        : requirment: matplotlib.pyplot
+        : requirement: matplotlib.pyplot
         """
         try:

src/sas/sascalc/calculator/resolution_calculator.py

@@ -1007 +1007 @@
         try:
             detector_offset = self.sample2detector_distance[1]
-        except Exception as ex:
-            logger.error(ex)
+        except:
+            logger.error(sys.exc_value)
 
         # detector size in [no of pix_x,no of pix_y]
@@ -1057 +1057 @@
         # qx_value and qy_value values in array
         qx_value = qx_value.repeat(detector_pix_nums_y)
-        qx_value = qx_value.reshape(int(detector_pix_nums_x), int(detector_pix_nums_y))
+        qx_value = qx_value.reshape(detector_pix_nums_x, detector_pix_nums_y)
         qy_value = qy_value.repeat(detector_pix_nums_x)
-        qy_value = qy_value.reshape(int(detector_pix_nums_y), int(detector_pix_nums_x))
+        qy_value = qy_value.reshape(detector_pix_nums_y, detector_pix_nums_x)
         qy_value = qy_value.transpose()
 
@@ -1094 +1094 @@
             output.qx_data = qx_value
             output.qy_data = qy_value
-        except Exception as ex:
-            logger.error(ex)
+        except:
+            logger.error(sys.exc_value)
 
         return output

src/sas/sascalc/calculator/sas_gen.py

@@ -118 +118 @@
         self.is_avg = is_avg
 
-    def _gen(self, x, y, i):
+    def _gen(self, qx, qy):
         """
         Evaluate the function
@@ -129 +129 @@
         pos_y = self.data_y
         pos_z = self.data_z
-        len_x = len(pos_x)
         if self.is_avg is None:
-            len_x *= -1
             pos_x, pos_y, pos_z = transform_center(pos_x, pos_y, pos_z)
-        len_q = len(x)
         sldn = copy.deepcopy(self.data_sldn)
         sldn -= self.params['solvent_SLD']
-        model = mod.new_GenI(len_x, pos_x, pos_y, pos_z,
-                             sldn, self.data_mx, self.data_my,
-                             self.data_mz, self.data_vol,
-                             self.params['Up_frac_in'],
-                             self.params['Up_frac_out'],
-                             self.params['Up_theta'])
-        if y == []:
-            mod.genicom(model, len_q, x, i)
-        else:
-            mod.genicomXY(model, len_q, x, y, i)
+        # **** WARNING **** new_GenI holds pointers to numpy vectors
+        # be sure that they are contiguous double precision arrays and make 
+        # sure the GC doesn't eat them before genicom is called.
+        # TODO: rewrite so that the parameters are passed directly to genicom
+        args = (
+            (1 if self.is_avg else 0),
+            pos_x, pos_y, pos_z,
+            sldn, self.data_mx, self.data_my,
+            self.data_mz, self.data_vol,
+            self.params['Up_frac_in'],
+            self.params['Up_frac_out'],
+            self.params['Up_theta'])
+        model = mod.new_GenI(*args)
+        if len(qy):
+            qx, qy = _vec(qx), _vec(qy)
+            I_out = np.empty_like(qx)
+            #print("npoints", qx.shape, "npixels", pos_x.shape)
+            mod.genicomXY(model, qx, qy, I_out)
+            #print("I_out after", I_out)
+        else:
+            qx = _vec(qx)
+            I_out = np.empty_like(qx)
+            mod.genicom(model, qx, I_out)
         vol_correction = self.data_total_volume / self.params['total_volume']
-        return  self.params['scale'] * vol_correction * i + \
-                        self.params['background']
+        result = (self.params['scale'] * vol_correction * I_out
+                  + self.params['background'])
+        return result
 
     def set_sld_data(self, sld_data=None):
@@ -156 +167 @@
         self.sld_data = sld_data
         self.data_pos_unit = sld_data.pos_unit
-        self.data_x = sld_data.pos_x
-        self.data_y = sld_data.pos_y
-        self.data_z = sld_data.pos_z
-        self.data_sldn = sld_data.sld_n
-        self.data_mx = sld_data.sld_mx
-        self.data_my = sld_data.sld_my
-        self.data_mz = sld_data.sld_mz
-        self.data_vol = sld_data.vol_pix
+        self.data_x = _vec(sld_data.pos_x)
+        self.data_y = _vec(sld_data.pos_y)
+        self.data_z = _vec(sld_data.pos_z)
+        self.data_sldn = _vec(sld_data.sld_n)
+        self.data_mx = _vec(sld_data.sld_mx)
+        self.data_my = _vec(sld_data.sld_my)
+        self.data_mz = _vec(sld_data.sld_mz)
+        self.data_vol = _vec(sld_data.vol_pix)
         self.data_total_volume = sum(sld_data.vol_pix)
         self.params['total_volume'] = sum(sld_data.vol_pix)
@@ -180 +191 @@
         :return: (I value)
         """
-        if x.__class__.__name__ == 'list':
+        if isinstance(x, list):
             if len(x[1]) > 0:
                 msg = "Not a 1D."
                 raise ValueError(msg)
-            i_out = np.zeros_like(x[0])
             # 1D I is found at y =0 in the 2D pattern
-            out = self._gen(x[0], [], i_out)
+            out = self._gen(x[0], [])
             return out
         else:
@@ -199 +209 @@
         :Use this runXY() for the computation
         """
-        if x.__class__.__name__ == 'list':
-            i_out = np.zeros_like(x[0])
-            out = self._gen(x[0], x[1], i_out)
-            return out
+        if isinstance(x, list):
+            return self._gen(x[0], x[1])
         else:
             msg = "Q must be given as list of qx's and qy's"
@@ -214 +222 @@
                       where qx,qy are 1D ndarrays (for 2D).
         """
-        if qdist.__class__.__name__ == 'list':
-            if len(qdist[1]) < 1:
-                out = self.run(qdist)
-            else:
-                out = self.runXY(qdist)
-            return out
+        if isinstance(qdist, list):
+            return self.run(qdist) if len(qdist[1]) < 1 else self.runXY(qdist)
         else:
             mesg = "evalDistribution is expecting an ndarray of "
             mesg += "a list [qx,qy] where qx,qy are arrays."
             raise RuntimeError(mesg)
+
+def _vec(v):
+    return np.ascontiguousarray(v, 'd')
 
 class OMF2SLD(object):
@@ -1041 +1048 @@
         self.line_z = line_z
 
+def _get_data_path(*path_parts):
+    from os.path import realpath, join as joinpath, dirname, abspath
+    # in sas/sascalc/calculator;  want sas/sasview/test
+    return joinpath(dirname(realpath(__file__)),
+                    '..', '..', 'sasview', 'test', *path_parts)
+
 def test_load():
     """
@@ -1046 +1059 @@
     """
     from mpl_toolkits.mplot3d import Axes3D
-    current_dir = os.path.abspath(os.path.curdir)
-    print(current_dir)
-    for i in range(6):
-        current_dir, _ = os.path.split(current_dir)
-        tfile = os.path.join(current_dir, "test", "CoreXY_ShellZ.txt")
-        ofile = os.path.join(current_dir, "test", "A_Raw_Example-1.omf")
-        if os.path.isfile(tfile):
-            tfpath = tfile
-            ofpath = ofile
-            break
+    tfpath = _get_data_path("1d_data", "CoreXY_ShellZ.txt")
+    ofpath = _get_data_path("coordinate_data", "A_Raw_Example-1.omf")
+    if not os.path.isfile(tfpath) or not os.path.isfile(ofpath):
+        raise ValueError("file(s) not found: %r, %r"%(tfpath, ofpath))
     reader = SLDReader()
     oreader = OMFReader()
-    output = decode(reader.read(tfpath))
-    ooutput = decode(oreader.read(ofpath))
+    output = reader.read(tfpath)
+    ooutput = oreader.read(ofpath)
     foutput = OMF2SLD()
     foutput.set_data(ooutput)
@@ -1088 +1095 @@
     plt.show()
 
+def test_save():
+    ofpath = _get_data_path("coordinate_data", "A_Raw_Example-1.omf")
+    if not os.path.isfile(ofpath):
+        raise ValueError("file(s) not found: %r"%(ofpath,))
+    oreader = OMFReader()
+    omfdata = oreader.read(ofpath)
+    omf2sld = OMF2SLD()
+    omf2sld.set_data(omfdata)
+    writer = SLDReader()
+    writer.write("out.txt", omf2sld.output)
+
 def test():
     """
         Test code
     """
-    current_dir = os.path.abspath(os.path.curdir)
-    for i in range(3):
-        current_dir, _ = os.path.split(current_dir)
-        ofile = os.path.join(current_dir, "test", "A_Raw_Example-1.omf")
-        if os.path.isfile(ofile):
-            ofpath = ofile
-            break
+    ofpath = _get_data_path("coordinate_data", "A_Raw_Example-1.omf")
+    if not os.path.isfile(ofpath):
+        raise ValueError("file(s) not found: %r"%(ofpath,))
     oreader = OMFReader()
-    ooutput = decode(oreader.read(ofpath))
-    foutput = OMF2SLD()
-    foutput.set_data(ooutput)
-    writer = SLDReader()
-    writer.write(os.path.join(os.path.dirname(ofpath), "out.txt"),
-                 foutput.output)
+    omfdata = oreader.read(ofpath)
+    omf2sld = OMF2SLD()
+    omf2sld.set_data(omfdata)
     model = GenSAS()
-    model.set_sld_data(foutput.output)
-    x = np.arange(1000)/10000. + 1e-5
-    y = np.arange(1000)/10000. + 1e-5
-    i = np.zeros(1000)
-    model.runXY([x, y, i])
+    model.set_sld_data(omf2sld.output)
+    x = np.linspace(0, 0.1, 11)[1:]
+    return model.runXY([x, x])
 
 if __name__ == "__main__":
+    #test_load()
+    #test_save()
+    #print(test())
     test()
-    test_load()