Changeset 9a23253e in sasview for pr_inversion/c_extensions/Cinvertor.c
- Timestamp:
- Jun 30, 2008 2:18:00 PM (16 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:
- a4bd2ac
- Parents:
- 7bb61da
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
pr_inversion/c_extensions/Cinvertor.c
r896abb3 r9a23253e 3 3 * Cinvertor is the base class for the Invertor class 4 4 * and provides the underlying computations. 5 * 5 * 6 6 */ 7 7 #include <Python.h> … … 38 38 * Cinvertor is the base class for the Invertor class 39 39 * and provides the underlying computations. 40 * 40 * 41 41 */ 42 42 typedef struct { 43 PyObject_HEAD 43 PyObject_HEAD 44 44 /// Internal data structure 45 Invertor_params params; 45 Invertor_params params; 46 46 } Cinvertor; 47 47 … … 51 51 { 52 52 invertor_dealloc(&(self->params)); 53 53 54 54 self->ob_type->tp_free((PyObject*)self); 55 55 … … 60 60 { 61 61 Cinvertor *self; 62 62 63 63 self = (Cinvertor *)type->tp_alloc(type, 0); 64 64 65 65 return (PyObject *)self; 66 66 } … … 69 69 Cinvertor_init(Cinvertor *self, PyObject *args, PyObject *kwds) 70 70 { 71 if (self != NULL) { 71 if (self != NULL) { 72 72 // Create parameters 73 73 invertor_init(&(self->params)); … … 82 82 }; 83 83 84 const char set_x_doc[] = 84 const char set_x_doc[] = 85 85 "Function to set the x data\n" 86 86 "Takes an array of doubles as input.\n" … … 97 97 double *data; 98 98 int i; 99 99 100 100 if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 101 101 OUTVECTOR(data_obj,data,ndata); 102 102 103 103 free(self->params.x); 104 104 self->params.x = (double*) malloc(ndata*sizeof(double)); 105 105 106 106 if(self->params.x==NULL) { 107 PyErr_SetString(CinvertorError, 107 PyErr_SetString(CinvertorError, 108 108 "Cinvertor.set_x: problem allocating memory."); 109 return NULL; 110 } 111 109 return NULL; 110 } 111 112 112 for (i=0; i<ndata; i++) { 113 113 self->params.x[i] = data[i]; 114 114 } 115 115 116 116 //self->params.x = data; 117 117 self->params.npoints = ndata; 118 return Py_BuildValue("i", self->params.npoints); 119 } 120 121 const char get_x_doc[] = 118 return Py_BuildValue("i", self->params.npoints); 119 } 120 121 const char get_x_doc[] = 122 122 "Function to get the x data\n" 123 123 "Takes an array of doubles as input.\n" … … 129 129 double *data; 130 130 int i; 131 131 132 132 if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 133 133 OUTVECTOR(data_obj, data, ndata); 134 134 135 135 // Check that the input array is large enough 136 136 if (ndata < self->params.npoints) { 137 PyErr_SetString(CinvertorError, 137 PyErr_SetString(CinvertorError, 138 138 "Cinvertor.get_x: input array too short for data."); 139 return NULL; 140 } 141 139 return NULL; 140 } 141 142 142 for(i=0; i<self->params.npoints; i++){ 143 143 data[i] = self->params.x[i]; 144 144 } 145 146 return Py_BuildValue("i", self->params.npoints); 147 } 148 149 const char set_y_doc[] = 145 146 return Py_BuildValue("i", self->params.npoints); 147 } 148 149 const char set_y_doc[] = 150 150 "Function to set the y data\n" 151 151 "Takes an array of doubles as input.\n" … … 162 162 double *data; 163 163 int i; 164 164 165 165 if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 166 166 OUTVECTOR(data_obj,data,ndata); 167 167 168 168 free(self->params.y); 169 169 self->params.y = (double*) malloc(ndata*sizeof(double)); 170 170 171 171 if(self->params.y==NULL) { 172 PyErr_SetString(CinvertorError, 172 PyErr_SetString(CinvertorError, 173 173 "Cinvertor.set_y: problem allocating memory."); 174 return NULL; 175 } 176 174 return NULL; 175 } 176 177 177 for (i=0; i<ndata; i++) { 178 178 self->params.y[i] = data[i]; 179 } 180 179 } 180 181 181 //self->params.y = data; 182 182 self->params.ny = ndata; 183 return Py_BuildValue("i", self->params.ny); 184 } 185 186 const char get_y_doc[] = 183 return Py_BuildValue("i", self->params.ny); 184 } 185 186 const char get_y_doc[] = 187 187 "Function to get the y data\n" 188 188 "Takes an array of doubles as input.\n" … … 194 194 double *data; 195 195 int i; 196 196 197 197 if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 198 198 OUTVECTOR(data_obj, data, ndata); 199 199 200 200 // Check that the input array is large enough 201 201 if (ndata < self->params.ny) { 202 PyErr_SetString(CinvertorError, 202 PyErr_SetString(CinvertorError, 203 203 "Cinvertor.get_y: input array too short for data."); 204 return NULL; 205 } 206 204 return NULL; 205 } 206 207 207 for(i=0; i<self->params.ny; i++){ 208 208 data[i] = self->params.y[i]; 209 209 } 210 211 return Py_BuildValue("i", self->params.npoints); 212 } 213 214 const char set_err_doc[] = 210 211 return Py_BuildValue("i", self->params.npoints); 212 } 213 214 const char set_err_doc[] = 215 215 "Function to set the err data\n" 216 216 "Takes an array of doubles as input.\n" … … 227 227 double *data; 228 228 int i; 229 229 230 230 if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 231 231 OUTVECTOR(data_obj,data,ndata); 232 232 233 233 free(self->params.err); 234 234 self->params.err = (double*) malloc(ndata*sizeof(double)); 235 235 236 236 if(self->params.err==NULL) { 237 PyErr_SetString(CinvertorError, 237 PyErr_SetString(CinvertorError, 238 238 "Cinvertor.set_err: problem allocating memory."); 239 return NULL; 240 } 241 239 return NULL; 240 } 241 242 242 for (i=0; i<ndata; i++) { 243 243 self->params.err[i] = data[i]; 244 244 } 245 245 246 246 //self->params.err = data; 247 247 self->params.nerr = ndata; 248 return Py_BuildValue("i", self->params.nerr); 249 } 250 251 const char get_err_doc[] = 248 return Py_BuildValue("i", self->params.nerr); 249 } 250 251 const char get_err_doc[] = 252 252 "Function to get the err data\n" 253 253 "Takes an array of doubles as input.\n" … … 259 259 double *data; 260 260 int i; 261 261 262 262 if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 263 263 OUTVECTOR(data_obj, data, ndata); 264 264 265 265 // Check that the input array is large enough 266 266 if (ndata < self->params.nerr) { 267 PyErr_SetString(CinvertorError, 267 PyErr_SetString(CinvertorError, 268 268 "Cinvertor.get_err: input array too short for data."); 269 return NULL; 270 } 271 269 return NULL; 270 } 271 272 272 for(i=0; i<self->params.nerr; i++){ 273 273 data[i] = self->params.err[i]; 274 274 } 275 276 return Py_BuildValue("i", self->params.npoints); 277 } 278 279 const char is_valid_doc[] = 275 276 return Py_BuildValue("i", self->params.npoints); 277 } 278 279 const char is_valid_doc[] = 280 280 "Check the validity of the stored data\n" 281 281 " @return: Returns the number of points if it's all good, -1 otherwise"; … … 291 291 } else { 292 292 return Py_BuildValue("i", -1); 293 } 294 } 295 296 const char set_dmax_doc[] = 293 } 294 } 295 296 const char set_has_bck_doc[] = 297 "Sets background flag\n"; 298 299 /** 300 * Sets the maximum distance 301 */ 302 static PyObject * set_has_bck(Cinvertor *self, PyObject *args) { 303 int has_bck; 304 305 if (!PyArg_ParseTuple(args, "i", &has_bck)) return NULL; 306 self->params.has_bck = has_bck; 307 return Py_BuildValue("i", self->params.has_bck); 308 } 309 310 const char get_has_bck_doc[] = 311 "Gets background flag\n"; 312 313 /** 314 * Gets the maximum distance 315 */ 316 static PyObject * get_has_bck(Cinvertor *self, PyObject *args) { 317 return Py_BuildValue("i", self->params.has_bck); 318 } 319 320 const char set_dmax_doc[] = 297 321 "Sets the maximum distance\n"; 298 322 … … 302 326 static PyObject * set_dmax(Cinvertor *self, PyObject *args) { 303 327 double d_max; 304 328 305 329 if (!PyArg_ParseTuple(args, "d", &d_max)) return NULL; 306 330 self->params.d_max = d_max; 307 return Py_BuildValue("d", self->params.d_max); 308 } 309 310 const char get_dmax_doc[] = 331 return Py_BuildValue("d", self->params.d_max); 332 } 333 334 const char get_dmax_doc[] = 311 335 "Gets the maximum distance\n"; 312 336 … … 315 339 */ 316 340 static PyObject * get_dmax(Cinvertor *self, PyObject *args) { 317 return Py_BuildValue("d", self->params.d_max); 318 } 319 320 const char set_qmin_doc[] = 341 return Py_BuildValue("d", self->params.d_max); 342 } 343 344 const char set_slit_height_doc[] = 345 "Sets the slit height in units of q [A-1]\n"; 346 347 /** 348 * Sets the slit height 349 */ 350 static PyObject * set_slit_height(Cinvertor *self, PyObject *args) { 351 double slit_height; 352 353 if (!PyArg_ParseTuple(args, "d", &slit_height)) return NULL; 354 self->params.slit_height = slit_height; 355 return Py_BuildValue("d", self->params.slit_height); 356 } 357 358 const char get_slit_height_doc[] = 359 "Gets the slit height\n"; 360 361 /** 362 * Gets the slit height 363 */ 364 static PyObject * get_slit_height(Cinvertor *self, PyObject *args) { 365 return Py_BuildValue("d", self->params.slit_height); 366 } 367 368 const char set_slit_width_doc[] = 369 "Sets the slit width in units of q [A-1]\n"; 370 371 /** 372 * Sets the slit width 373 */ 374 static PyObject * set_slit_width(Cinvertor *self, PyObject *args) { 375 double slit_width; 376 377 if (!PyArg_ParseTuple(args, "d", &slit_width)) return NULL; 378 self->params.slit_width = slit_width; 379 return Py_BuildValue("d", self->params.slit_width); 380 } 381 382 const char get_slit_width_doc[] = 383 "Gets the slit width\n"; 384 385 /** 386 * Gets the slit width 387 */ 388 static PyObject * get_slit_width(Cinvertor *self, PyObject *args) { 389 return Py_BuildValue("d", self->params.slit_width); 390 } 391 392 393 const char set_qmin_doc[] = 321 394 "Sets the minimum q\n"; 322 395 … … 326 399 static PyObject * set_qmin(Cinvertor *self, PyObject *args) { 327 400 double q_min; 328 401 329 402 if (!PyArg_ParseTuple(args, "d", &q_min)) return NULL; 330 403 self->params.q_min = q_min; 331 return Py_BuildValue("d", self->params.q_min); 332 } 333 334 const char get_qmin_doc[] = 404 return Py_BuildValue("d", self->params.q_min); 405 } 406 407 const char get_qmin_doc[] = 335 408 "Gets the minimum q\n"; 336 409 … … 339 412 */ 340 413 static PyObject * get_qmin(Cinvertor *self, PyObject *args) { 341 return Py_BuildValue("d", self->params.q_min); 342 } 343 344 const char set_qmax_doc[] = 414 return Py_BuildValue("d", self->params.q_min); 415 } 416 417 const char set_qmax_doc[] = 345 418 "Sets the maximum q\n"; 346 419 … … 350 423 static PyObject * set_qmax(Cinvertor *self, PyObject *args) { 351 424 double q_max; 352 425 353 426 if (!PyArg_ParseTuple(args, "d", &q_max)) return NULL; 354 427 self->params.q_max = q_max; 355 return Py_BuildValue("d", self->params.q_max); 356 } 357 358 const char get_qmax_doc[] = 428 return Py_BuildValue("d", self->params.q_max); 429 } 430 431 const char get_qmax_doc[] = 359 432 "Gets the maximum q\n"; 360 433 … … 363 436 */ 364 437 static PyObject * get_qmax(Cinvertor *self, PyObject *args) { 365 return Py_BuildValue("d", self->params.q_max); 366 } 367 368 const char set_alpha_doc[] = 438 return Py_BuildValue("d", self->params.q_max); 439 } 440 441 const char set_alpha_doc[] = 369 442 "Sets the alpha parameter\n"; 370 443 371 444 static PyObject * set_alpha(Cinvertor *self, PyObject *args) { 372 445 double alpha; 373 446 374 447 if (!PyArg_ParseTuple(args, "d", &alpha)) return NULL; 375 448 self->params.alpha = alpha; 376 return Py_BuildValue("d", self->params.alpha); 377 } 378 379 const char get_alpha_doc[] = 449 return Py_BuildValue("d", self->params.alpha); 450 } 451 452 const char get_alpha_doc[] = 380 453 "Gets the alpha parameter\n"; 381 454 … … 384 457 */ 385 458 static PyObject * get_alpha(Cinvertor *self, PyObject *args) { 386 return Py_BuildValue("d", self->params.alpha); 387 } 388 389 const char get_nx_doc[] = 459 return Py_BuildValue("d", self->params.alpha); 460 } 461 462 const char get_nx_doc[] = 390 463 "Gets the number of x points\n"; 391 464 … … 394 467 */ 395 468 static PyObject * get_nx(Cinvertor *self, PyObject *args) { 396 return Py_BuildValue("i", self->params.npoints); 397 } 398 399 const char get_ny_doc[] = 469 return Py_BuildValue("i", self->params.npoints); 470 } 471 472 const char get_ny_doc[] = 400 473 "Gets the number of y points\n"; 401 474 … … 404 477 */ 405 478 static PyObject * get_ny(Cinvertor *self, PyObject *args) { 406 return Py_BuildValue("i", self->params.ny); 407 } 408 409 const char get_nerr_doc[] = 479 return Py_BuildValue("i", self->params.ny); 480 } 481 482 const char get_nerr_doc[] = 410 483 "Gets the number of err points\n"; 411 484 … … 414 487 */ 415 488 static PyObject * get_nerr(Cinvertor *self, PyObject *args) { 416 return Py_BuildValue("i", self->params.nerr); 417 } 418 419 420 const char residuals_doc[] = 489 return Py_BuildValue("i", self->params.nerr); 490 } 491 492 493 const char residuals_doc[] = 421 494 "Function to call to evaluate the residuals\n" 422 495 "for P(r) inversion\n" … … 441 514 // Number of slices in regularization term estimate 442 515 int nslice = 25; 443 516 444 517 PyObject *data_obj; 445 518 Py_ssize_t npars; 446 447 if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 448 519 520 if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 521 449 522 OUTVECTOR(data_obj,pars,npars); 450 523 451 524 // PyList of residuals 452 525 // Should create this list only once and refill it … … 454 527 455 528 regterm = reg_term(pars, self->params.d_max, npars, nslice); 456 529 457 530 for(i=0; i<self->params.npoints; i++) { 458 531 diff = self->params.y[i] - iq(pars, self->params.d_max, npars, self->params.x[i]); 459 532 residual = diff*diff / (self->params.err[i]*self->params.err[i]); 460 533 tmp = residual; 461 534 462 535 // regularization term 463 536 residual += self->params.alpha * regterm; 464 537 465 538 if (PyList_SetItem(residuals, i, Py_BuildValue("d",residual) ) < 0){ 466 PyErr_SetString(CinvertorError, 539 PyErr_SetString(CinvertorError, 467 540 "Cinvertor.residuals: error setting residual."); 468 541 return NULL; 469 542 }; 470 543 471 544 } 472 545 473 546 return residuals; 474 547 } 475 548 476 const char pr_residuals_doc[] = 549 const char pr_residuals_doc[] = 477 550 "Function to call to evaluate the residuals\n" 478 551 "for P(r) minimization (for testing purposes)\n" … … 498 571 // Number of slices in regularization term estimate 499 572 int nslice = 25; 500 573 501 574 PyObject *data_obj; 502 575 Py_ssize_t npars; 503 504 if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 505 576 577 if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 578 506 579 OUTVECTOR(data_obj,pars,npars); 507 580 508 581 // Should create this list only once and refill it 509 582 residuals = PyList_New(self->params.npoints); … … 511 584 regterm = reg_term(pars, self->params.d_max, npars, nslice); 512 585 513 586 514 587 for(i=0; i<self->params.npoints; i++) { 515 588 diff = self->params.y[i] - pr(pars, self->params.d_max, npars, self->params.x[i]); 516 589 residual = diff*diff / (self->params.err[i]*self->params.err[i]); 517 590 tmp = residual; 518 591 519 592 // regularization term 520 593 residual += self->params.alpha * regterm; 521 594 522 595 if (PyList_SetItem(residuals, i, Py_BuildValue("d",residual) ) < 0){ 523 PyErr_SetString(CinvertorError, 596 PyErr_SetString(CinvertorError, 524 597 "Cinvertor.residuals: error setting residual."); 525 598 return NULL; 526 599 }; 527 600 528 601 } 529 602 530 603 return residuals; 531 604 } 532 605 533 const char get_iq_doc[] = 606 const char get_iq_doc[] = 534 607 "Function to call to evaluate the scattering intensity\n" 535 608 " @param args: c-parameters, and q\n" … … 546 619 PyObject *data_obj; 547 620 Py_ssize_t npars; 548 621 549 622 if (!PyArg_ParseTuple(args, "Od", &data_obj, &q)) return NULL; 550 623 OUTVECTOR(data_obj,pars,npars); 551 624 552 625 iq_value = iq(pars, self->params.d_max, npars, q); 553 return Py_BuildValue("f", iq_value); 554 } 555 556 const char get_pr_doc[] = 626 return Py_BuildValue("f", iq_value); 627 } 628 629 const char get_pr_doc[] = 557 630 "Function to call to evaluate P(r)\n" 558 631 " @param args: c-parameters and r\n" … … 569 642 PyObject *data_obj; 570 643 Py_ssize_t npars; 571 644 572 645 if (!PyArg_ParseTuple(args, "Od", &data_obj, &r)) return NULL; 573 646 OUTVECTOR(data_obj,pars,npars); 574 647 575 648 pr_value = pr(pars, self->params.d_max, npars, r); 576 return Py_BuildValue("f", pr_value); 577 } 578 579 const char get_pr_err_doc[] = 649 return Py_BuildValue("f", pr_value); 650 } 651 652 const char get_pr_err_doc[] = 580 653 "Function to call to evaluate P(r) with errors\n" 581 654 " @param args: c-parameters and r\n" … … 596 669 PyObject *err_obj; 597 670 Py_ssize_t npars2; 598 int i; 599 671 int i; 672 600 673 if (!PyArg_ParseTuple(args, "OOd", &data_obj, &err_obj, &r)) return NULL; 601 OUTVECTOR(data_obj,pars,npars); 674 OUTVECTOR(data_obj,pars,npars); 602 675 OUTVECTOR(err_obj,pars_err,npars2); 603 676 604 677 pr_err(pars, pars_err, self->params.d_max, npars, r, &pr_value, &pr_err_value); 605 return Py_BuildValue("ff", pr_value, pr_err_value); 606 } 607 608 const char basefunc_ft_doc[] = 678 return Py_BuildValue("ff", pr_value, pr_err_value); 679 } 680 681 const char basefunc_ft_doc[] = 609 682 "Returns the value of the nth Fourier transofrmed base function\n" 610 683 " @param args: c-parameters, n and q\n" … … 614 687 double d_max, q; 615 688 int n; 616 689 617 690 if (!PyArg_ParseTuple(args, "did", &d_max, &n, &q)) return NULL; 618 return Py_BuildValue("f", ortho_transformed(d_max, n, q)); 619 620 } 621 622 const char oscillations_doc[] = 691 return Py_BuildValue("f", ortho_transformed(d_max, n, q)); 692 693 } 694 695 const char oscillations_doc[] = 623 696 "Returns the value of the oscillation figure of merit for\n" 624 697 "the given set of coefficients. For a sphere, the oscillation\n" … … 632 705 Py_ssize_t npars; 633 706 double oscill, norm; 634 707 635 708 if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 636 709 OUTVECTOR(data_obj,pars,npars); 637 710 638 711 oscill = reg_term(pars, self->params.d_max, npars, 100); 639 712 norm = int_p2(pars, self->params.d_max, npars, 100); 640 return Py_BuildValue("f", sqrt(oscill/norm)/acos(-1.0)*self->params.d_max ); 641 642 } 643 644 const char get_peaks_doc[] = 713 return Py_BuildValue("f", sqrt(oscill/norm)/acos(-1.0)*self->params.d_max ); 714 715 } 716 717 const char get_peaks_doc[] = 645 718 "Returns the number of peaks in the output P(r) distrubution\n" 646 719 "for the given set of coefficients.\n" … … 653 726 Py_ssize_t npars; 654 727 int count; 655 728 656 729 if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 657 730 OUTVECTOR(data_obj,pars,npars); 658 731 659 732 count = npeaks(pars, self->params.d_max, npars, 100); 660 733 661 return Py_BuildValue("i", count ); 662 663 } 664 665 const char get_positive_doc[] = 734 return Py_BuildValue("i", count ); 735 736 } 737 738 const char get_positive_doc[] = 666 739 "Returns the fraction of P(r) that is positive over\n" 667 740 "the full range of r for the given set of coefficients.\n" … … 674 747 Py_ssize_t npars; 675 748 double fraction; 676 749 677 750 if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 678 751 OUTVECTOR(data_obj,pars,npars); 679 752 680 753 fraction = positive_integral(pars, self->params.d_max, npars, 100); 681 754 682 return Py_BuildValue("f", fraction ); 683 684 } 685 686 const char get_pos_err_doc[] = 755 return Py_BuildValue("f", fraction ); 756 757 } 758 759 const char get_pos_err_doc[] = 687 760 "Returns the fraction of P(r) that is 1 standard deviation\n" 688 761 "above zero over the full range of r for the given set of coefficients.\n" … … 698 771 Py_ssize_t npars2; 699 772 double fraction; 700 773 701 774 if (!PyArg_ParseTuple(args, "OO", &data_obj, &err_obj)) return NULL; 702 OUTVECTOR(data_obj,pars,npars); 775 OUTVECTOR(data_obj,pars,npars); 703 776 OUTVECTOR(err_obj,pars_err,npars2); 704 777 705 778 fraction = positive_errors(pars, pars_err, self->params.d_max, npars, 51); 706 779 707 return Py_BuildValue("f", fraction ); 708 709 } 710 780 return Py_BuildValue("f", fraction ); 781 782 } 783 784 const char get_rg_doc[] = 785 "Returns the value of the radius of gyration Rg.\n" 786 " @param args: c-parameters\n" 787 " @return: Rg"; 788 789 static PyObject * get_rg(Cinvertor *self, PyObject *args) { 790 double *pars; 791 double *pars_err; 792 PyObject *data_obj; 793 Py_ssize_t npars; 794 Py_ssize_t npars2; 795 double value; 796 797 if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 798 OUTVECTOR(data_obj,pars,npars); 799 800 value = rg(pars, self->params.d_max, npars, 101); 801 802 return Py_BuildValue("f", value ); 803 804 } 805 806 const char get_iq0_doc[] = 807 "Returns the value of I(q=0).\n" 808 " @param args: c-parameters\n" 809 " @return: I(q=0)"; 810 811 static PyObject * get_iq0(Cinvertor *self, PyObject *args) { 812 double *pars; 813 double *pars_err; 814 PyObject *data_obj; 815 Py_ssize_t npars; 816 Py_ssize_t npars2; 817 double value; 818 819 if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 820 OUTVECTOR(data_obj,pars,npars); 821 822 value = 4.0*acos(-1.0)*int_pr(pars, self->params.d_max, npars, 101); 823 824 return Py_BuildValue("f", value ); 825 826 } 827 828 /** 829 * Check whether a q-value is within acceptabel limits 830 * Return 1 if accepted, 0 if rejected. 831 */ 832 int accept_q(Cinvertor *self, double q) { 833 if (self->params.q_min>0 && q<self->params.q_min) return 0; 834 if (self->params.q_max>0 && q>self->params.q_max) return 0; 835 return 1; 836 } 837 838 const char get_matrix_doc[] = 839 "Returns A matrix and b vector for least square problem.\n" 840 " @param nfunc: number of base functions\n" 841 " @param nr: number of r-points used when evaluating reg term.\n" 842 " @param a: A array to fill\n" 843 " @param b: b vector to fill\n" 844 " @return: 0"; 845 846 static PyObject * get_matrix(Cinvertor *self, PyObject *args) { 847 double *a; 848 double *b; 849 PyObject *a_obj; 850 PyObject *b_obj; 851 Py_ssize_t n_a; 852 Py_ssize_t n_b; 853 // Number of bins for regularization term evaluation 854 int nr, nfunc; 855 int i, j, i_r; 856 double r, sqrt_alpha, pi; 857 double tmp; 858 int offset; 859 860 if (!PyArg_ParseTuple(args, "iiOO", &nfunc, &nr, &a_obj, &b_obj)) return NULL; 861 OUTVECTOR(a_obj,a,n_a); 862 OUTVECTOR(b_obj,b,n_b); 863 864 assert(n_b>=nfunc); 865 assert(n_a>=nfunc*(nr*self->params.npoints)); 866 867 sqrt_alpha = sqrt(self->params.alpha); 868 pi = acos(-1.0); 869 offset = (self->params.has_bck==1) ? 0 : 1; 870 871 for (j=0; j<nfunc; j++) { 872 for (i=0; i<self->params.npoints; i++) { 873 if (accept_q(self, self->params.x[i])){ 874 if (self->params.has_bck==1 && j==0) { 875 a[i*nfunc+j] = 1.0/self->params.err[i]; 876 } else { 877 a[i*nfunc+j] = ortho_transformed(self->params.d_max, j+offset, self->params.x[i])/self->params.err[i]; 878 } 879 } 880 } 881 for (i_r=0; i_r<nr; i_r++){ 882 if (self->params.has_bck==1 && j==0) { 883 a[(i_r+self->params.npoints)*nfunc+j] = 0.0; 884 } else { 885 r = self->params.d_max/nr*i_r; 886 tmp = pi*(j+offset)/self->params.d_max; 887 a[(i_r+self->params.npoints)*nfunc+j] = sqrt_alpha * 1.0/nr*self->params.d_max*2.0* 888 (2.0*pi*(j+offset)/self->params.d_max*cos(pi*(j+offset)*r/self->params.d_max) + 889 tmp*tmp*r * sin(pi*(j+offset)*r/self->params.d_max)); 890 } 891 } 892 } 893 894 for (i=0; i<self->params.npoints; i++) { 895 if (accept_q(self, self->params.x[i])){ 896 b[i] = self->params.y[i]/self->params.err[i]; 897 } 898 } 899 900 return Py_BuildValue("i", 0); 901 902 } 903 904 const char get_invcov_matrix_doc[] = 905 " Compute the inverse covariance matrix, defined as inv_cov = a_transposed x a.\n" 906 " @param nfunc: number of base functions\n" 907 " @param nr: number of r-points used when evaluating reg term.\n" 908 " @param a: A array to fill\n" 909 " @param inv_cov: inverse covariance array to be filled\n" 910 " @return: 0"; 911 912 static PyObject * get_invcov_matrix(Cinvertor *self, PyObject *args) { 913 double *a; 914 PyObject *a_obj; 915 Py_ssize_t n_a; 916 double *inv_cov; 917 PyObject *cov_obj; 918 Py_ssize_t n_cov; 919 int nr, nfunc; 920 int i, j, k; 921 922 if (!PyArg_ParseTuple(args, "iiOO", &nfunc, &nr, &a_obj, &cov_obj)) return NULL; 923 OUTVECTOR(a_obj,a,n_a); 924 OUTVECTOR(cov_obj,inv_cov,n_cov); 925 926 assert(n_cov>=nfunc*nfunc); 927 assert(n_a>=nfunc*(nr+self->params.npoints)); 928 929 for (i=0; i<nfunc; i++) { 930 for (j=0; j<nfunc; j++) { 931 inv_cov[i*nfunc+j] = 0.0; 932 for (k=0; k<nr+self->params.npoints; k++) { 933 inv_cov[i*nfunc+j] += a[k*nfunc+i]*a[k*nfunc+j]; 934 } 935 } 936 } 937 return Py_BuildValue("i", 0); 938 } 939 940 const char get_reg_size_doc[] = 941 " Compute the covariance matrix, defined as inv_cov = a_transposed x a.\n" 942 " @param nfunc: number of base functions\n" 943 " @param nr: number of r-points used when evaluating reg term.\n" 944 " @param a: A array to fill\n" 945 " @param inv_cov: inverse covariance array to be filled\n" 946 " @return: 0"; 947 948 static PyObject * get_reg_size(Cinvertor *self, PyObject *args) { 949 double *a; 950 PyObject *a_obj; 951 Py_ssize_t n_a; 952 int nr, nfunc; 953 int i, j; 954 double sum_sig, sum_reg; 955 956 if (!PyArg_ParseTuple(args, "iiO", &nfunc, &nr, &a_obj)) return NULL; 957 OUTVECTOR(a_obj,a,n_a); 958 959 assert(n_a>=nfunc*(nr+self->params.npoints)); 960 961 sum_sig = 0.0; 962 sum_reg = 0.0; 963 for (j=0; j<nfunc; j++){ 964 for (i=0; i<self->params.npoints; i++){ 965 if (accept_q(self, self->params.x[i])==1) 966 sum_sig += (a[i*nfunc+j])*(a[i*nfunc+j]); 967 } 968 for (i=0; i<nr; i++){ 969 sum_reg += (a[(i+self->params.npoints)*nfunc+j])*(a[(i+self->params.npoints)*nfunc+j]); 970 } 971 } 972 return Py_BuildValue("ff", sum_sig, sum_reg); 973 } 711 974 712 975 const char eeeget_qmin_doc[] = "\ … … 714 977 \n\ 715 978 This is the second line."; 716 const char eeeset_qmin_doc[] = 979 const char eeeset_qmin_doc[] = 717 980 "This is a multiline doc string.\n" 718 981 "\n" … … 736 999 {"set_alpha", (PyCFunction)set_alpha, METH_VARARGS, set_alpha_doc}, 737 1000 {"get_alpha", (PyCFunction)get_alpha, METH_VARARGS, get_alpha_doc}, 1001 {"set_slit_width", (PyCFunction)set_slit_width, METH_VARARGS, set_slit_width_doc}, 1002 {"get_slit_width", (PyCFunction)get_slit_width, METH_VARARGS, get_slit_width_doc}, 1003 {"set_slit_height", (PyCFunction)set_slit_height, METH_VARARGS, set_slit_height_doc}, 1004 {"get_slit_height", (PyCFunction)get_slit_height, METH_VARARGS, get_slit_height_doc}, 1005 {"set_has_bck", (PyCFunction)set_has_bck, METH_VARARGS, set_has_bck_doc}, 1006 {"get_has_bck", (PyCFunction)get_has_bck, METH_VARARGS, get_has_bck_doc}, 738 1007 {"get_nx", (PyCFunction)get_nx, METH_VARARGS, get_nx_doc}, 739 1008 {"get_ny", (PyCFunction)get_ny, METH_VARARGS, get_ny_doc}, … … 748 1017 {"get_positive", (PyCFunction)get_positive, METH_VARARGS, get_positive_doc}, 749 1018 {"get_pos_err", (PyCFunction)get_pos_err, METH_VARARGS, get_pos_err_doc}, 750 1019 {"rg", (PyCFunction)get_rg, METH_VARARGS, get_rg_doc}, 1020 {"iq0", (PyCFunction)get_iq0, METH_VARARGS, get_iq0_doc}, 1021 {"_get_matrix", (PyCFunction)get_matrix, METH_VARARGS, get_matrix_doc}, 1022 {"_get_invcov_matrix", (PyCFunction)get_invcov_matrix, METH_VARARGS, get_invcov_matrix_doc}, 1023 {"_get_reg_size", (PyCFunction)get_reg_size, METH_VARARGS, get_reg_size_doc}, 1024 751 1025 {NULL} 752 1026 }; … … 796 1070 797 1071 static PyMethodDef module_methods[] = { 798 {NULL} 1072 {NULL} 799 1073 }; 800 1074 … … 802 1076 * Function used to add the model class to a module 803 1077 * @param module: module to add the class to 804 */ 1078 */ 805 1079 void addCinvertor(PyObject *module) { 806 1080 PyObject *d; 807 1081 808 1082 if (PyType_Ready(&CinvertorType) < 0) 809 1083 return; … … 811 1085 Py_INCREF(&CinvertorType); 812 1086 PyModule_AddObject(module, "Cinvertor", (PyObject *)&CinvertorType); 813 1087 814 1088 d = PyModule_GetDict(module); 815 1089 CinvertorError = PyErr_NewException("Cinvertor.error", NULL, NULL); … … 822 1096 #endif 823 1097 PyMODINIT_FUNC 824 initpr_inversion(void) 1098 initpr_inversion(void) 825 1099 { 826 1100 PyObject* m; … … 828 1102 m = Py_InitModule3("pr_inversion", module_methods, 829 1103 "C extension module for inversion to P(r)."); 830 1104 831 1105 addCinvertor(m); 832 1106 }
Note: See TracChangeset
for help on using the changeset viewer.