Changeset 54b0650 in sasview


Ignore:
Timestamp:
Nov 7, 2017 1:05:12 PM (7 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, magnetic_scatt, release-4.2.2, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
Children:
f926abb
Parents:
0957bb3a
Message:

fix C interface to sldi after py3 conversion

Location:
src/sas/sascalc/calculator
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r0957bb3a r54b0650  
    2525 * @param s_theta: angle (from x-axis) of the up spin in degree 
    2626 */ 
    27 void initGenI(GenI* this, int npix, double* x, double* y, double* z, double* sldn, 
     27void initGenI(GenI* this, int is_avg, int npix, double* x, double* y, double* z, double* sldn, 
    2828                        double* mx, double* my, double* mz, double* voli, 
    2929                        double in_spin, double out_spin, 
    3030                        double s_theta) { 
     31        this->is_avg = is_avg; 
    3132        this->n_pix = npix; 
    3233        this->x_val = x; 
    … …  
    7374        // Loop over q-values and multiply apply matrix 
    7475 
     76        //printf("npoints: %d, npix: %d\n", npoints, this->n_pix); 
    7577        for(int i=0; i<npoints; i++){ 
    7678                //I_out[i] = 0.0; 
    … …  
    150152        // Assumes that q doesn't have qz component and sld_n is all real 
    151153        //double Pi = 4.0*atan(1.0); 
    152         int is_sym = this->n_pix < 0; 
    153154        double qr = 0.0; 
    154155        double sumj; 
    155156        double sld_j = 0.0; 
    156157        double count = 0.0; 
    157         int n_pix = is_sym ? -this->n_pix : this->n_pix; 
    158158        //Assume that pixel volumes are given in vol_pix in A^3 unit 
    159159        // Loop over q-values and multiply apply matrix 
    160160        for(int i=0; i<npoints; i++){ 
    161161                sumj =0.0; 
    162                 for(int j=0; j<n_pix; j++){ 
     162                for(int j=0; j<this->n_pix; j++){ 
    163163                        //Isotropic: Assumes all slds are real (no magnetic) 
    164164                        //Also assumes there is no polarization: No dependency on spin 
    165                         if (is_sym == 1){ 
     165                        if (this->is_avg == 1){ 
    166166                                // approximation for a spherical symmetric particle 
    167167                                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]; 
    … …  
    177177                                //full calculation 
    178178                                //pragma omp parallel for 
    179                                 for(int k=0; k<n_pix; k++){ 
     179                                for(int k=0; k<this->n_pix; k++){ 
    180180                                        sld_j =  this->sldn_val[j] * this->sldn_val[k] * this->vol_pix[j] * this->vol_pix[k]; 
    181181                                        qr = (this->x_val[j]-this->x_val[k])*(this->x_val[j]-this->x_val[k])+ 
    … …  
    196196                } 
    197197                I_out[i] = sumj; 
    198                 if (is_sym == 1){ 
     198                if (this->is_avg == 1) { 
    199199                        I_out[i] *= sumj; 
    200200                } 

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

    r0957bb3a r54b0650  
    1010typedef struct { 
    1111        // vectors 
     12        int is_avg; 
    1213        int n_pix; 
    1314        double* x_val; 
    … …  
    2627 
    2728// Constructor 
    28 void initGenI(GenI*, int npix, double* x, double* y, double* z, double* sldn, 
    29                 double* mx, double* my, double* mz, double* voli, 
     29void initGenI(GenI*, int is_avg, int npix, double* x, double* y, double* z, 
     30                double* sldn, double* mx, double* my, double* mz, double* voli, 
    3031                double in_spin, double out_spin, 
    3132                double s_theta); 

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

    r0957bb3a r54b0650  
    3535void 
    3636del_sld2i(PyObject *obj){ 
     37#if PY_MAJOR_VERSION < 3 
     38        GenI* sld2i = (GenI *)obj; 
     39#else 
    3740        GenI* sld2i = (GenI *)(PyCapsule_GetPointer(obj, "GenI")); 
     41#endif 
    3842        PyMem_Free((void *)sld2i); 
    3943} 
    … …  
    5155        PyObject *mz_val_obj; 
    5256        PyObject *vol_pix_obj; 
    53         Py_ssize_t n_x; 
    54         //PyObject rlimit_obj; 
    55         //PyObject npoints_obj; 
    56         //PyObject nrbins_obj; 
    57         //PyObject nphibins_obj; 
    58         int n_pix; 
     57        Py_ssize_t n_x, n_y, n_z, n_sld, n_mx, n_my, n_mz, n_vol_pix; 
     58        int is_avg; 
    5959        double* x_val; 
    6060        double* y_val; 
    … …  
    6969        double stheta; 
    7070 
    71         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; 
    72         OUTVECTOR(x_val_obj, x_val, n_x); 
    73         OUTVECTOR(y_val_obj, y_val, n_x); 
    74         OUTVECTOR(z_val_obj, z_val, n_x); 
    75         OUTVECTOR(sldn_val_obj, sldn_val, n_x); 
    76         OUTVECTOR(mx_val_obj, mx_val, n_x); 
    77         OUTVECTOR(my_val_obj, my_val, n_x); 
    78         OUTVECTOR(mz_val_obj, mz_val, n_x); 
    79         OUTVECTOR(vol_pix_obj, vol_pix, n_x); 
    80         GenI* sld2i =  PyMem_Malloc(sizeof(GenI)); 
     71        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; 
     72        INVECTOR(x_val_obj, x_val, n_x); 
     73        INVECTOR(y_val_obj, y_val, n_y); 
     74        INVECTOR(z_val_obj, z_val, n_z); 
     75        INVECTOR(sldn_val_obj, sldn_val, n_sld); 
     76        INVECTOR(mx_val_obj, mx_val, n_mx); 
     77        INVECTOR(my_val_obj, my_val, n_my); 
     78        INVECTOR(mz_val_obj, mz_val, n_mz); 
     79        INVECTOR(vol_pix_obj, vol_pix, n_vol_pix); 
     80        GenI* sld2i = PyMem_Malloc(sizeof(GenI)); 
     81        //printf("sldi:%p\n", sld2i); 
    8182        if (sld2i != NULL) { 
    82                 initGenI(sld2i, n_pix,x_val,y_val,z_val,sldn_val,mx_val,my_val,mz_val,vol_pix,inspin,outspin,stheta); 
     83                initGenI(sld2i,is_avg,n_x,x_val,y_val,z_val,sldn_val,mx_val,my_val,mz_val,vol_pix,inspin,outspin,stheta); 
    8384        } 
    84         return PyCapsule_New(sld2i, "GenI", del_sld2i); 
     85        PyObject *obj = PyCapsule_New(sld2i, "GenI", del_sld2i); 
     86        //printf("constructed %p\n", obj); 
     87        return obj; 
    8588} 
    8689 
    … …  
    8992 */ 
    9093PyObject * genicom_inputXY(PyObject *self, PyObject *args) { 
    91         int npoints; 
     94        PyObject *gen_obj; 
    9295        PyObject *qx_obj; 
     96        PyObject *qy_obj; 
     97        PyObject *I_out_obj; 
     98        Py_ssize_t n_qx, n_qy, n_out; 
    9399        double *qx; 
    94         PyObject *qy_obj; 
    95100        double *qy; 
    96         PyObject *I_out_obj; 
    97         Py_ssize_t n_out; 
    98101        double *I_out; 
    99         PyObject *gen_obj; 
    100102 
    101         if (!PyArg_ParseTuple(args, "OiOOO",  &gen_obj, &npoints, &qx_obj, &qy_obj, &I_out_obj)) return NULL; 
    102         OUTVECTOR(qx_obj, qx, n_out); 
    103         OUTVECTOR(qy_obj, qy, n_out); 
     103        if (!PyArg_ParseTuple(args, "OOOO",  &gen_obj, &qx_obj, &qy_obj, &I_out_obj)) return NULL; 
     104        GenI* sld2i = (GenI *)PyCapsule_GetPointer(gen_obj, "GenI"); 
     105        INVECTOR(qx_obj, qx, n_qx); 
     106        INVECTOR(qy_obj, qy, n_qy); 
    104107        OUTVECTOR(I_out_obj, I_out, n_out); 
    105108 
    106109        // Sanity check 
    107         //if(n_in!=n_out) return Py_BuildValue("i",-1); 
     110        //if(n_q!=n_out) return Py_BuildValue("i",-1); 
    108111 
    109         // Set the array pointers 
    110         GenI* sld2i = (GenI *)PyCapsule_GetPointer(gen_obj, "GenI"); 
    111  
    112         genicomXY(sld2i, npoints, qx, qy, I_out); 
     112        genicomXY(sld2i, n_qx, qx, qy, I_out); 
    113113        //return PyCObject_FromVoidPtr(s, del_genicom); 
    114114        return Py_BuildValue("i",1); 
    … …  
    119119 */ 
    120120PyObject * genicom_input(PyObject *self, PyObject *args) { 
    121         int npoints; 
     121        PyObject *gen_obj; 
    122122        PyObject *q_obj; 
     123        PyObject *I_out_obj; 
     124        Py_ssize_t n_q, n_out; 
    123125        double *q; 
    124         PyObject *I_out_obj; 
    125         Py_ssize_t n_out; 
    126126        double *I_out; 
    127         PyObject *gen_obj; 
    128127 
    129         if (!PyArg_ParseTuple(args, "OiOO",  &gen_obj, &npoints, &q_obj, &I_out_obj)) return NULL; 
    130         OUTVECTOR(q_obj, q, n_out); 
     128        if (!PyArg_ParseTuple(args, "OOO",  &gen_obj, &q_obj, &I_out_obj)) return NULL; 
     129        GenI *sld2i = (GenI *)PyCapsule_GetPointer(gen_obj, "GenI"); 
     130        INVECTOR(q_obj, q, n_q); 
    131131        OUTVECTOR(I_out_obj, I_out, n_out); 
    132132 
    133133        // Sanity check 
    134         //if(n_in!=n_out) return Py_BuildValue("i",-1); 
     134        //if (n_q!=n_out) return Py_BuildValue("i",-1); 
    135135 
    136         // Set the array pointers 
    137         GenI *sld2i = (GenI *)PyCapsule_GetPointer(gen_obj, "GenI"); 
    138  
    139         genicom(sld2i, npoints, q, I_out); 
    140         //return PyCObject_FromVoidPtr(s, del_genicom); 
     136        genicom(sld2i, n_q, q, I_out); 
    141137        return Py_BuildValue("i",1); 
    142138} 

  • src/sas/sascalc/calculator/sas_gen.py

    rb58265c3 r54b0650  
    118118        self.is_avg = is_avg 
    119119 
    120     def _gen(self, x, y, i): 
     120    def _gen(self, qx, qy): 
    121121        """ 
    122122        Evaluate the function 
    … …  
    129129        pos_y = self.data_y 
    130130        pos_z = self.data_z 
    131         len_x = len(pos_x) 
    132131        if self.is_avg is None: 
    133             len_x *= -1 
    134132            pos_x, pos_y, pos_z = transform_center(pos_x, pos_y, pos_z) 
    135         len_q = len(x) 
    136133        sldn = copy.deepcopy(self.data_sldn) 
    137134        sldn -= self.params['solvent_SLD'] 
    138         model = mod.new_GenI(len_x, pos_x, pos_y, pos_z, 
     135        model = mod.new_GenI((1 if self.is_avg else 0), 
     136                             pos_x, pos_y, pos_z, 
    139137                             sldn, self.data_mx, self.data_my, 
    140138                             self.data_mz, self.data_vol, 
    … …  
    142140                             self.params['Up_frac_out'], 
    143141                             self.params['Up_theta']) 
    144         if y == []: 
    145             mod.genicom(model, len_q, x, i) 
    146         else: 
    147             mod.genicomXY(model, len_q, x, y, i) 
     142        if len(qy): 
     143            qx, qy = _vec(qx), _vec(qy) 
     144            I_out = np.empty_like(qx) 
     145            mod.genicomXY(model, qx, qy, I_out) 
     146        else: 
     147            qx = _vec(qx) 
     148            I_out = np.empty_like(qx) 
     149            mod.genicom(model, qx, I_out) 
    148150        vol_correction = self.data_total_volume / self.params['total_volume'] 
    149         return  self.params['scale'] * vol_correction * i + \ 
    150                         self.params['background'] 
     151        result = (self.params['scale'] * vol_correction * I_out 
     152                  + self.params['background']) 
     153        return result 
    151154 
    152155    def set_sld_data(self, sld_data=None): 
    … …  
    156159        self.sld_data = sld_data 
    157160        self.data_pos_unit = sld_data.pos_unit 
    158         self.data_x = sld_data.pos_x 
    159         self.data_y = sld_data.pos_y 
    160         self.data_z = sld_data.pos_z 
    161         self.data_sldn = sld_data.sld_n 
    162         self.data_mx = sld_data.sld_mx 
    163         self.data_my = sld_data.sld_my 
    164         self.data_mz = sld_data.sld_mz 
    165         self.data_vol = sld_data.vol_pix 
     161        self.data_x = _vec(sld_data.pos_x) 
     162        self.data_y = _vec(sld_data.pos_y) 
     163        self.data_z = _vec(sld_data.pos_z) 
     164        self.data_sldn = _vec(sld_data.sld_n) 
     165        self.data_mx = _vec(sld_data.sld_mx) 
     166        self.data_my = _vec(sld_data.sld_my) 
     167        self.data_mz = _vec(sld_data.sld_mz) 
     168        self.data_vol = _vec(sld_data.vol_pix) 
    166169        self.data_total_volume = sum(sld_data.vol_pix) 
    167170        self.params['total_volume'] = sum(sld_data.vol_pix) 
    … …  
    180183        :return: (I value) 
    181184        """ 
    182         if x.__class__.__name__ == 'list': 
     185        if isinstance(x, list): 
    183186            if len(x[1]) > 0: 
    184187                msg = "Not a 1D." 
    185188                raise ValueError(msg) 
    186             i_out = np.zeros_like(x[0]) 
    187189            # 1D I is found at y =0 in the 2D pattern 
    188             out = self._gen(x[0], [], i_out) 
     190            out = self._gen(x[0], []) 
    189191            return out 
    190192        else: 
    … …  
    199201        :Use this runXY() for the computation 
    200202        """ 
    201         if x.__class__.__name__ == 'list': 
    202             i_out = np.zeros_like(x[0]) 
    203             out = self._gen(x[0], x[1], i_out) 
    204             return out 
     203        if isinstance(x, list): 
     204            return self._gen(x[0], x[1]) 
    205205        else: 
    206206            msg = "Q must be given as list of qx's and qy's" 
    … …  
    214214                      where qx,qy are 1D ndarrays (for 2D). 
    215215        """ 
    216         if qdist.__class__.__name__ == 'list': 
    217             if len(qdist[1]) < 1: 
    218                 out = self.run(qdist) 
    219             else: 
    220                 out = self.runXY(qdist) 
    221             return out 
     216        if isinstance(qdist, list): 
     217            return self.run(qdist) if len(qdist[1]) < 1 else self.runXY(qdist) 
    222218        else: 
    223219            mesg = "evalDistribution is expecting an ndarray of " 
    224220            mesg += "a list [qx,qy] where qx,qy are arrays." 
    225221            raise RuntimeError(mesg) 
     222 
     223def _vec(v): 
     224    return np.ascontiguousarray(v, 'd') 
    226225 
    227226class OMF2SLD(object): 
    … …  
    10411040        self.line_z = line_z 
    10421041 
     1042def _get_data_path(*path_parts): 
     1043    from os.path import realpath, join as joinpath, dirname, abspath 
     1044    # in sas/sascalc/calculator;  want sas/sasview/test 
     1045    return joinpath(dirname(realpath(__file__)), 
     1046                    '..', '..', 'sasview', 'test', *path_parts) 
     1047 
    10431048def test_load(): 
    10441049    """ 
    … …  
    10461051    """ 
    10471052    from mpl_toolkits.mplot3d import Axes3D 
    1048     current_dir = os.path.abspath(os.path.curdir) 
    1049     print(current_dir) 
    1050     for i in range(6): 
    1051         current_dir, _ = os.path.split(current_dir) 
    1052         tfile = os.path.join(current_dir, "test", "CoreXY_ShellZ.txt") 
    1053         ofile = os.path.join(current_dir, "test", "A_Raw_Example-1.omf") 
    1054         if os.path.isfile(tfile): 
    1055             tfpath = tfile 
    1056             ofpath = ofile 
    1057             break 
     1053    tfpath = _get_data_path("1d_data", "CoreXY_ShellZ.txt") 
     1054    ofpath = _get_data_path("coordinate_data", "A_Raw_Example-1.omf") 
     1055    if not os.path.isfile(tfpath) or not os.path.isfile(ofpath): 
     1056        raise ValueError("file(s) not found: %r, %r"%(tfpath, ofpath)) 
    10581057    reader = SLDReader() 
    10591058    oreader = OMFReader() 
    … …  
    10921091        Test code 
    10931092    """ 
    1094     current_dir = os.path.abspath(os.path.curdir) 
    1095     for i in range(3): 
    1096         current_dir, _ = os.path.split(current_dir) 
    1097         ofile = os.path.join(current_dir, "test", "A_Raw_Example-1.omf") 
    1098         if os.path.isfile(ofile): 
    1099             ofpath = ofile 
    1100             break 
     1093    ofpath = _get_data_path("coordinate_data", "A_Raw_Example-1.omf") 
     1094    if not os.path.isfile(ofpath): 
     1095        raise ValueError("file(s) not found: %r"%(ofpath,)) 
    11011096    oreader = OMFReader() 
    11021097    ooutput = decode(oreader.read(ofpath)) 
    … …  
    11041099    foutput.set_data(ooutput) 
    11051100    writer = SLDReader() 
    1106     writer.write(os.path.join(os.path.dirname(ofpath), "out.txt"), 
    1107                  foutput.output) 
     1101    writer.write("out.txt", foutput.output) 
    11081102    model = GenSAS() 
    11091103    model.set_sld_data(foutput.output) 
    … …  
    11141108 
    11151109if __name__ == "__main__": 
     1110    #test_load() 
    11161111    test() 
    1117     test_load() 
Note: See TracChangeset for help on using the changeset viewer.