Changeset 9bd69098 in sasview for sansmodels/src/sans/models/c_models/CEllipticalCylinderModel.cpp
- Timestamp:
- Aug 7, 2009 9:10:34 AM (15 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:
- 58c6ba6
- Parents:
- 83a25da
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sansmodels/src/sans/models/c_models/CEllipticalCylinderModel.cpp
r0f5bc9f r9bd69098 22 22 * 23 23 */ 24 #define NO_IMPORT_ARRAY 25 #define PY_ARRAY_UNIQUE_SYMBOL PyArray_API_sans 24 26 25 27 extern "C" { 26 28 #include <Python.h> 29 #include <arrayobject.h> 27 30 #include "structmember.h" 28 31 #include <stdio.h> … … 150 153 } 151 154 } 152 153 154 155 /** 155 156 * Function to call to evaluate model 156 * @param args: input q or [q,phi]157 * @return: function value157 * @param args: input numpy array q[] 158 * @return: numpy array object 158 159 */ 159 static PyObject * run(CEllipticalCylinderModel *self, PyObject *args) { 160 double q_value, phi_value; 161 PyObject* pars; 162 int npars; 160 161 static PyObject *evaluateOneDim(EllipticalCylinderModel* model, PyArrayObject *q){ 162 PyArrayObject *result; 163 164 // Check validity of array q , q must be of dimension 1, an array of double 165 if (q->nd != 1 || q->descr->type_num != PyArray_DOUBLE) 166 { 167 //const char * message= "Invalid array: q->nd=%d,type_num=%d\n",q->nd,q->descr->type_num; 168 //PyErr_SetString(PyExc_ValueError , message); 169 return NULL; 170 } 171 result = (PyArrayObject *)PyArray_FromDims(q->nd, (int *)(q->dimensions), 172 PyArray_DOUBLE); 173 if (result == NULL) { 174 const char * message= "Could not create result "; 175 PyErr_SetString(PyExc_RuntimeError , message); 176 return NULL; 177 } 178 for (int i = 0; i < q->dimensions[0]; i++){ 179 double q_value = *(double *)(q->data + i*q->strides[0]); 180 double *result_value = (double *)(result->data + i*result->strides[0]); 181 *result_value =(*model)(q_value); 182 } 183 return PyArray_Return(result); 184 } 185 /** 186 * Function to call to evaluate model 187 * @param args: input numpy array [q[],phi[]] 188 * @return: numpy array object 189 */ 190 static PyObject * evaluateTwoDim( EllipticalCylinderModel* model, 191 PyArrayObject *q, PyArrayObject *phi) 192 { 193 PyArrayObject *result; 194 //check validity of input vectors 195 if (q->nd != 1 || q->descr->type_num != PyArray_DOUBLE 196 || phi->nd != 1 || phi->descr->type_num != PyArray_DOUBLE 197 || phi->dimensions[0] != q->dimensions[0]){ 198 199 //const char * message= "Invalid array: q->nd=%d,type_num=%d\n",q->nd,q->descr->type_num; 200 PyErr_SetString(PyExc_ValueError ,"wrong input"); 201 return NULL; 202 } 203 result= (PyArrayObject *)PyArray_FromDims(q->nd,(int*)(q->dimensions), PyArray_DOUBLE); 204 205 if (result == NULL){ 206 const char * message= "Could not create result "; 207 PyErr_SetString(PyExc_RuntimeError , message); 208 return NULL; 209 } 210 211 for (int i = 0; i < q->dimensions[0]; i++) { 212 double q_value = *(double *)(q->data + i*q->strides[0]); 213 double phi_value = *(double *)(phi->data + i*phi->strides[0]); 214 double *result_value = (double *)(result->data + i*result->strides[0]); 215 if (q_value == 0) 216 *result_value = 0.0; 217 else 218 *result_value = model->evaluate_rphi(q_value, phi_value); 219 } 220 return PyArray_Return(result); 221 } 222 /** 223 * Function to call to evaluate model 224 * @param args: input numpy array [x[],y[]] 225 * @return: numpy array object 226 */ 227 static PyObject * evaluateTwoDimXY( EllipticalCylinderModel* model, 228 PyArrayObject *x, PyArrayObject *y) 229 { 230 PyArrayObject *result; 231 int i,j, x_len, y_len, dims[2]; 232 //check validity of input vectors 233 if (x->nd != 2 || x->descr->type_num != PyArray_DOUBLE 234 || y->nd != 2 || y->descr->type_num != PyArray_DOUBLE 235 || y->dimensions[1] != x->dimensions[0]){ 236 const char * message= "evaluateTwoDimXY expect 2 numpy arrays"; 237 PyErr_SetString(PyExc_ValueError , message); 238 return NULL; 239 } 240 241 if (PyArray_Check(x) && PyArray_Check(y)) { 242 x_len = dims[0]= x->dimensions[0]; 243 y_len = dims[1]= y->dimensions[1]; 244 245 // Make a new double matrix of same dims 246 result=(PyArrayObject *) PyArray_FromDims(2,dims,NPY_DOUBLE); 247 if (result == NULL){ 248 const char * message= "Could not create result "; 249 PyErr_SetString(PyExc_RuntimeError , message); 250 return NULL; 251 } 252 253 /* Do the calculation. */ 254 for ( i=0; i< x_len; i++) { 255 for ( j=0; j< y_len; j++) { 256 double x_value = *(double *)(x->data + i*x->strides[0]); 257 double y_value = *(double *)(y->data + j*y->strides[1]); 258 double *result_value = (double *)(result->data + 259 i*result->strides[0] + j*result->strides[1]); 260 *result_value = (*model)(x_value, y_value); 261 } 262 } 263 return PyArray_Return(result); 264 265 }else{ 266 PyErr_SetString(CEllipticalCylinderModelError, 267 "CEllipticalCylinderModel.evaluateTwoDimXY couldn't run."); 268 return NULL; 269 } 270 } 271 /** 272 * evalDistribution function evaluate a model function with input vector 273 * @param args: input q as vector or [qx, qy] where qx, qy are vectors 274 * 275 */ 276 static PyObject * evalDistribution(CEllipticalCylinderModel *self, PyObject *args){ 277 PyObject *qx, *qy; 278 PyArrayObject * pars; 279 int npars ,mpars; 163 280 164 281 // Get parameters … … 194 311 if ( !PyArg_ParseTuple(args,"O",&pars) ) { 195 312 PyErr_SetString(CEllipticalCylinderModelError, 196 "CEllipticalCylinderModel. run expects a q value.");313 "CEllipticalCylinderModel.evalDistribution expects a q value."); 197 314 return NULL; 198 315 } 199 200 // Check params201 if( PyList_Check(pars)==1) {316 // Check params 317 318 if(PyArray_Check(pars)==1) { 202 319 203 // Length of list should be 2 for I(q,phi) 204 npars = PyList_GET_SIZE(pars); 205 if(npars!=2) { 320 // Length of list should 1 or 2 321 npars = pars->nd; 322 if(npars==1) { 323 // input is a numpy array 324 if (PyArray_Check(pars)) { 325 return evaluateOneDim(self->model, (PyArrayObject*)pars); 326 } 327 }else{ 328 PyErr_SetString(CEllipticalCylinderModelError, 329 "CEllipticalCylinderModel.evalDistribution expect numpy array of one dimension."); 330 return NULL; 331 } 332 }else if( PyList_Check(pars)==1) { 333 // Length of list should be 2 for I(qx,qy) 334 mpars = PyList_GET_SIZE(pars); 335 if(mpars!=2) { 206 336 PyErr_SetString(CEllipticalCylinderModelError, 207 "CEllipticalCylinderModel. run expects a double ora list of dimension 2.");337 "CEllipticalCylinderModel.evalDistribution expects a list of dimension 2."); 208 338 return NULL; 209 339 } 210 // We have a vector q, get the q and phi values at which 211 // to evaluate I(q,phi) 212 q_value = CEllipticalCylinderModel_readDouble(PyList_GET_ITEM(pars,0)); 213 phi_value = CEllipticalCylinderModel_readDouble(PyList_GET_ITEM(pars,1)); 214 // Skip zero 215 if (q_value==0) { 216 return Py_BuildValue("d",0.0); 217 } 218 return Py_BuildValue("d",(*(self->model)).evaluate_rphi(q_value,phi_value)); 219 220 } else { 221 222 // We have a scalar q, we will evaluate I(q) 223 q_value = CEllipticalCylinderModel_readDouble(pars); 224 225 return Py_BuildValue("d",(*(self->model))(q_value)); 226 } 340 qx = PyList_GET_ITEM(pars,0); 341 qy = PyList_GET_ITEM(pars,1); 342 if (PyArray_Check(qx) && PyArray_Check(qy)) { 343 return evaluateTwoDimXY(self->model, (PyArrayObject*)qx, 344 (PyArrayObject*)qy); 345 }else{ 346 PyErr_SetString(CEllipticalCylinderModelError, 347 "CEllipticalCylinderModel.evalDistribution expect 2 numpy arrays in list."); 348 return NULL; 349 } 350 }else{ 351 PyErr_SetString(CEllipticalCylinderModelError, 352 "CEllipticalCylinderModel.evalDistribution couln't be run."); 353 return NULL; 354 } 227 355 } 228 356 229 357 /** 230 * Function to call to evaluate model in cartesian coordinates231 * @param args: input q or [q x, qy]]358 * Function to call to evaluate model 359 * @param args: input q or [q,phi] 232 360 * @return: function value 233 361 */ 234 static PyObject * run XY(CEllipticalCylinderModel *self, PyObject *args) {235 double q x_value, qy_value;362 static PyObject * run(CEllipticalCylinderModel *self, PyObject *args) { 363 double q_value, phi_value; 236 364 PyObject* pars; 237 365 int npars; … … 276 404 if( PyList_Check(pars)==1) { 277 405 406 // Length of list should be 2 for I(q,phi) 407 npars = PyList_GET_SIZE(pars); 408 if(npars!=2) { 409 PyErr_SetString(CEllipticalCylinderModelError, 410 "CEllipticalCylinderModel.run expects a double or a list of dimension 2."); 411 return NULL; 412 } 413 // We have a vector q, get the q and phi values at which 414 // to evaluate I(q,phi) 415 q_value = CEllipticalCylinderModel_readDouble(PyList_GET_ITEM(pars,0)); 416 phi_value = CEllipticalCylinderModel_readDouble(PyList_GET_ITEM(pars,1)); 417 // Skip zero 418 if (q_value==0) { 419 return Py_BuildValue("d",0.0); 420 } 421 return Py_BuildValue("d",(*(self->model)).evaluate_rphi(q_value,phi_value)); 422 423 } else { 424 425 // We have a scalar q, we will evaluate I(q) 426 q_value = CEllipticalCylinderModel_readDouble(pars); 427 428 return Py_BuildValue("d",(*(self->model))(q_value)); 429 } 430 } 431 432 /** 433 * Function to call to evaluate model in cartesian coordinates 434 * @param args: input q or [qx, qy]] 435 * @return: function value 436 */ 437 static PyObject * runXY(CEllipticalCylinderModel *self, PyObject *args) { 438 double qx_value, qy_value; 439 PyObject* pars; 440 int npars; 441 442 // Get parameters 443 444 // Reader parameter dictionary 445 self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") ); 446 self->model->cyl_psi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "cyl_psi") ); 447 self->model->length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "length") ); 448 self->model->r_minor = PyFloat_AsDouble( PyDict_GetItemString(self->params, "r_minor") ); 449 self->model->cyl_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "cyl_theta") ); 450 self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") ); 451 self->model->r_ratio = PyFloat_AsDouble( PyDict_GetItemString(self->params, "r_ratio") ); 452 self->model->contrast = PyFloat_AsDouble( PyDict_GetItemString(self->params, "contrast") ); 453 self->model->cyl_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "cyl_phi") ); 454 // Read in dispersion parameters 455 PyObject* disp_dict; 456 DispersionVisitor* visitor = new DispersionVisitor(); 457 disp_dict = PyDict_GetItemString(self->dispersion, "r_minor"); 458 self->model->r_minor.dispersion->accept_as_destination(visitor, self->model->r_minor.dispersion, disp_dict); 459 disp_dict = PyDict_GetItemString(self->dispersion, "r_ratio"); 460 self->model->r_ratio.dispersion->accept_as_destination(visitor, self->model->r_ratio.dispersion, disp_dict); 461 disp_dict = PyDict_GetItemString(self->dispersion, "length"); 462 self->model->length.dispersion->accept_as_destination(visitor, self->model->length.dispersion, disp_dict); 463 disp_dict = PyDict_GetItemString(self->dispersion, "cyl_theta"); 464 self->model->cyl_theta.dispersion->accept_as_destination(visitor, self->model->cyl_theta.dispersion, disp_dict); 465 disp_dict = PyDict_GetItemString(self->dispersion, "cyl_phi"); 466 self->model->cyl_phi.dispersion->accept_as_destination(visitor, self->model->cyl_phi.dispersion, disp_dict); 467 disp_dict = PyDict_GetItemString(self->dispersion, "cyl_psi"); 468 self->model->cyl_psi.dispersion->accept_as_destination(visitor, self->model->cyl_psi.dispersion, disp_dict); 469 470 471 // Get input and determine whether we have to supply a 1D or 2D return value. 472 if ( !PyArg_ParseTuple(args,"O",&pars) ) { 473 PyErr_SetString(CEllipticalCylinderModelError, 474 "CEllipticalCylinderModel.run expects a q value."); 475 return NULL; 476 } 477 478 // Check params 479 if( PyList_Check(pars)==1) { 480 278 481 // Length of list should be 2 for I(qx, qy)) 279 482 npars = PyList_GET_SIZE(pars); … … 350 553 {"runXY", (PyCFunction)runXY , METH_VARARGS, 351 554 "Evaluate the model at a given Q or Qx, Qy"}, 555 556 {"evalDistribution", (PyCFunction)evalDistribution , METH_VARARGS, 557 "Evaluate the model at a given Q or Qx, Qy vector "}, 352 558 {"reset", (PyCFunction)reset , METH_VARARGS, 353 559 "Reset pair correlation"}, … … 400 606 401 607 402 static PyMethodDef module_methods[] = {403 {NULL}404 };608 //static PyMethodDef module_methods[] = { 609 // {NULL} 610 //}; 405 611 406 612 /**
Note: See TracChangeset
for help on using the changeset viewer.