Changeset 9a23253e in sasview
- 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
- Location:
- pr_inversion
- Files:
-
- 6 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 } -
pr_inversion/c_extensions/invertor.c
rc61228f r9a23253e 20 20 pars->q_min = -1.0; 21 21 pars->q_max = -1.0; 22 pars->has_bck = 0; 22 23 } 23 24 … … 25 26 /** 26 27 * P(r) of a sphere, for test purposes 27 * 28 * 28 29 * @param R: radius of the sphere 29 30 * @param r: distance, in the same units as the radius … … 41 42 * Orthogonal functions: 42 43 * B(r) = 2r sin(pi*nr/d) 43 * 44 * 44 45 */ 45 46 double ortho(double d_max, int n, double r) { … … 49 50 /** 50 51 * Fourier transform of the nth orthogonal function 51 * 52 * 52 53 */ 53 54 double ortho_transformed(double d_max, int n, double q) { … … 57 58 58 59 /** 60 * Slit-smeared Fourier transform of the nth orthogonal function. 61 * Smearing follows Lake, Acta Cryst. (1967) 23, 191. 62 */ 63 double ortho_transformed_smeared(double d_max, int n, double height, double width, double q, int npts) { 64 double sum, value, y, z; 65 int i, j; 66 double fnpts; 67 sum = 0.0; 68 fnpts = (float)npts-1.0; 69 70 for(i=0; i<npts; i++) { 71 y = -width/2.0+width/fnpts*(float)i; 72 for(j=0; j<npts; j++) { 73 z = height/fnpts*(float)j; 74 sum += ortho_transformed(d_max, n, sqrt((q-y)*(q-y)+z*z)); 75 } 76 } 77 78 return sum/npts/npts/height/width; 79 } 80 81 /** 59 82 * First derivative in of the orthogonal function dB(r)/dr 60 * 83 * 61 84 */ 62 85 double ortho_derived(double d_max, int n, double r) { … … 80 103 */ 81 104 double pr(double *pars, double d_max, int n_c, double r) { 82 double sum = 0.0; 105 double sum = 0.0; 83 106 int i; 84 107 for (i=0; i<n_c; i++) { … … 91 114 * P(r) calculated from the expansion, with errors 92 115 */ 93 void pr_err(double *pars, double *err, double d_max, int n_c, 116 void pr_err(double *pars, double *err, double d_max, int n_c, 94 117 double r, double *pr_value, double *pr_value_err) { 95 double sum = 0.0; 118 double sum = 0.0; 96 119 double sum_err = 0.0; 97 120 double func_value; … … 109 132 *pr_value_err = sum; 110 133 } 111 } 134 } 112 135 113 136 /** … … 115 138 */ 116 139 double dprdr(double *pars, double d_max, int n_c, double r) { 117 double sum = 0.0; 118 int i; 140 double sum = 0.0; 141 int i; 119 142 for (i=0; i<n_c; i++) { 120 143 sum += pars[i] * 2.0*(sin(pi*(i+1)*r/d_max) + pi*(i+1)*r/d_max * cos(pi*(i+1)*r/d_max)); … … 127 150 */ 128 151 double reg_term(double *pars, double d_max, int n_c, int nslice) { 129 double sum = 0.0; 152 double sum = 0.0; 130 153 double r; 131 154 double deriv; … … 143 166 */ 144 167 double int_p2(double *pars, double d_max, int n_c, int nslice) { 145 double sum = 0.0; 146 double r; 168 double sum = 0.0; 169 double r; 147 170 double value; 148 171 int i; … … 156 179 157 180 /** 181 * Integral of P(r) 182 */ 183 double int_pr(double *pars, double d_max, int n_c, int nslice) { 184 double sum = 0.0; 185 double r; 186 double value; 187 int i; 188 for (i=0; i<nslice; i++) { 189 r = d_max/(1.0*nslice)*i; 190 value = pr(pars, d_max, n_c, r); 191 sum += value; 192 } 193 return sum/(1.0*nslice)*d_max; 194 } 195 196 /** 158 197 * Get the number of P(r) peaks. 159 198 */ 160 199 int npeaks(double *pars, double d_max, int n_c, int nslice) { 161 double r; 200 double r; 162 201 double value; 163 202 int i; … … 187 226 */ 188 227 double positive_integral(double *pars, double d_max, int n_c, int nslice) { 189 double r; 228 double r; 190 229 double value; 191 230 int i; 192 231 double sum_pos = 0.0; 193 232 double sum = 0.0; 194 233 195 234 for (i=0; i<nslice; i++) { 196 235 r = d_max/(1.0*nslice)*i; … … 207 246 */ 208 247 double positive_errors(double *pars, double *err, double d_max, int n_c, int nslice) { 209 double r; 210 double value; 211 int i; 248 double r; 249 double value; 250 int i; 212 251 double sum_pos = 0.0; 213 252 double sum = 0.0; 214 253 double pr_val; 215 254 double pr_val_err; 216 255 217 256 for (i=0; i<nslice; i++) { 218 257 r = d_max/(1.0*nslice)*i; … … 220 259 if (pr_val>pr_val_err) sum_pos += pr_val; 221 260 sum += fabs(pr_val); 222 261 223 262 224 263 } 225 264 return sum_pos/sum; 226 265 } 266 267 /** 268 * R_g radius of gyration calculation 269 * 270 * R_g**2 = integral[r**2 * p(r) dr] / (2.0 * integral[p(r) dr]) 271 */ 272 double rg(double *pars, double d_max, int n_c, int nslice) { 273 double sum_r2 = 0.0; 274 double sum = 0.0; 275 double r; 276 double value; 277 int i; 278 for (i=0; i<nslice; i++) { 279 r = d_max/(1.0*nslice)*i; 280 value = pr(pars, d_max, n_c, r); 281 sum += value; 282 sum_r2 += r*r*value; 283 } 284 return sqrt(sum_r2/(2.0*sum)); 285 } 286 -
pr_inversion/c_extensions/invertor.h
r896abb3 r9a23253e 15 15 double *err; 16 16 /// Number of q points 17 int npoints; 17 int npoints; 18 18 /// Number of I(q) points 19 int ny; 19 int ny; 20 20 /// Number of dI(q) points 21 int nerr; 21 int nerr; 22 22 /// Alpha value 23 23 double alpha; … … 26 26 /// Maximum q to include in inversion 27 27 double q_max; 28 } Invertor_params; 28 /// Flag for whether or not to evalute a constant background while inverting 29 int has_bck; 30 /// Slit height in units of q [A-1] 31 double slit_height; 32 /// Slit width in units of q [A-1] 33 double slit_width; 34 } Invertor_params; 29 35 30 36 void invertor_dealloc(Invertor_params *pars); 31 37 32 38 void invertor_init(Invertor_params *pars); 33 34 39 35 40 double pr_sphere(double R, double r); … … 42 47 double reg_term(double *pars, double d_max, int n_c, int nslice); 43 48 double int_p2(double *pars, double d_max, int n_c, int nslice); 44 void pr_err(double *pars, double *err, double d_max, int n_c, 49 void pr_err(double *pars, double *err, double d_max, int n_c, 45 50 double r, double *pr_value, double *pr_value_err); 46 51 int npeaks(double *pars, double d_max, int n_c, int nslice); 47 52 double positive_integral(double *pars, double d_max, int n_c, int nslice); 48 53 double positive_errors(double *pars, double *err, double d_max, int n_c, int nslice); 54 double rg(double *pars, double d_max, int n_c, int nslice); 55 double int_pr(double *pars, double d_max, int n_c, int nslice); 56 double ortho_transformed_smeared(double d_max, int n, double heigth, double width, double q, int npts); 49 57 50 58 #endif -
pr_inversion/invertor.py
rb00b487 r9a23253e 63 63 - get_pos_err(pars): returns the fraction of P(r) that is 1-sigma above zero 64 64 """ 65 #TODO: Allow for slit smearing. Smear each base function once before filling 66 # the A matrix. 67 65 68 ## Chisqr of the last computation 66 69 chi2 = 0 … … 75 78 ## Last errors on output values 76 79 cov = None 80 ## Flag to allow I(q) data with constant background 81 #has_bck = False 82 ## Background value 83 background = 0 84 77 85 78 86 def __init__(self): … … 106 114 elif name=='alpha': 107 115 return self.set_alpha(value) 116 elif name=='slit_height': 117 return self.set_slit_height(value) 118 elif name=='slit_width': 119 return self.set_slit_width(value) 120 elif name=='has_bck': 121 if value==True: 122 return self.set_has_bck(1) 123 elif value==False: 124 return self.set_has_bck(0) 125 else: 126 raise ValueError, "Invertor: has_bck can only be True or False" 108 127 109 128 return Cinvertor.__setattr__(self, name, value) … … 142 161 elif name=='alpha': 143 162 return self.get_alpha() 163 elif name=='slit_height': 164 return self.get_slit_height() 165 elif name=='slit_width': 166 return self.get_slit_width() 167 elif name=='has_bck': 168 value = self.get_has_bck() 169 if value==1: 170 return True 171 else: 172 return False 144 173 elif name in self.__dict__: 145 174 return self.__dict__[name] … … 161 190 invertor.y = self.y 162 191 invertor.err = self.err 192 invertor.has_bck = self.has_bck 163 193 164 194 return invertor … … 194 224 @return: c_out, c_cov - the coefficients with covariance matrix 195 225 """ 196 #TODO: call the pyhton implementation for now. In the future, translate this to C. 226 # Reset the background value before proceeding 227 self.background = 0.0 197 228 return self.lstsq(nfunc, nr=nr) 229 230 def iq(self, out, q): 231 """ 232 Function to call to evaluate the scattering intensity 233 @param args: c-parameters, and q 234 @return: I(q) 235 """ 236 return Cinvertor.iq(self, out, q)+self.background 198 237 199 238 def invert_optimize(self, nfunc=10, nr=20): … … 325 364 326 365 """ 327 import math 366 #TODO: Allow for background by starting at n=0 (since the base function 367 # is zero for n=0). 368 import math, time 328 369 from scipy.linalg.basic import lstsq 329 370 … … 339 380 if sqrt_alpha<0.0: 340 381 nq = 0 341 382 383 # If we need to fit the background, add a term 384 if self.has_bck==True: 385 nfunc_0 = nfunc 386 nfunc += 1 387 342 388 a = numpy.zeros([npts+nq, nfunc]) 343 389 b = numpy.zeros(npts+nq) 344 390 err = numpy.zeros([nfunc, nfunc]) 345 391 346 for j in range(nfunc): 347 for i in range(npts): 348 if self._accept_q(self.x[i]): 349 a[i][j] = self.basefunc_ft(self.d_max, j+1, self.x[i])/self.err[i] 350 351 #TODO: refactor this: i_q should really be i_r 352 for i_q in range(nq): 353 r = self.d_max/nq*i_q 354 #a[i_q+npts][j] = sqrt_alpha * 1.0/nq*self.d_max*2.0*math.fabs(math.sin(math.pi*(j+1)*r/self.d_max) + math.pi*(j+1)*r/self.d_max * math.cos(math.pi*(j+1)*r/self.d_max)) 355 a[i_q+npts][j] = sqrt_alpha * 1.0/nq*self.d_max*2.0*(2.0*math.pi*(j+1)/self.d_max*math.cos(math.pi*(j+1)*r/self.d_max) + math.pi**2*(j+1)**2*r/self.d_max**2 * math.sin(math.pi*(j+1)*r/self.d_max)) 356 357 for i in range(npts): 358 if self._accept_q(self.x[i]): 359 b[i] = self.y[i]/self.err[i] 360 392 # Construct the a matrix and b vector that represent the problem 393 self._get_matrix(nfunc, nq, a, b) 394 395 # Perform the inversion (least square fit) 361 396 c, chi2, rank, n = lstsq(a, b) 362 397 # Sanity check … … 367 402 self.chi2 = chi2 368 403 369 at = numpy.transpose(a)370 404 inv_cov = numpy.zeros([nfunc,nfunc]) 371 for i in range(nfunc): 372 for j in range(nfunc): 373 inv_cov[i][j] = 0.0 374 for k in range(npts+nr): 375 #if self._accept_q(self.x[i]): 376 inv_cov[i][j] += at[i][k]*a[k][j] 405 # Get the covariance matrix, defined as inv_cov = a_transposed * a 406 self._get_invcov_matrix(nfunc, nr, a, inv_cov) 377 407 378 408 # Compute the reg term size for the output 379 sum_sig = 0.0 380 sum_reg = 0.0 381 for j in range(nfunc): 382 for i in range(npts): 383 if self._accept_q(self.x[i]): 384 sum_sig += (a[i][j])**2 385 for i in range(nq): 386 sum_reg += (a[i+npts][j])**2 409 sum_sig, sum_reg = self._get_reg_size(nfunc, nr, a) 387 410 388 411 if math.fabs(self.alpha)>0: … … 403 426 404 427 # Keep a copy of the last output 405 self.out = c 406 self.cov = err 407 408 return c, err 428 if self.has_bck==False: 429 self.background = 0 430 self.out = c 431 self.cov = err 432 else: 433 self.background = c[0] 434 435 err_0 = numpy.zeros([nfunc, nfunc]) 436 c_0 = numpy.zeros(nfunc) 437 438 for i in range(nfunc_0): 439 c_0[i] = c[i+1] 440 for j in range(nfunc_0): 441 err_0[i][j] = err[i+1][j+1] 442 443 self.out = c_0 444 self.cov = err_0 445 446 return self.out, self.cov 447 448 def lstsq_bck(self, nfunc=5, nr=20): 449 #TODO: Allow for background by starting at n=0 (since the base function 450 # is zero for n=0). 451 import math 452 from scipy.linalg.basic import lstsq 453 454 if self.is_valid()<0: 455 raise RuntimeError, "Invertor: invalid data; incompatible data lengths." 456 457 self.nfunc = nfunc 458 # a -- An M x N matrix. 459 # b -- An M x nrhs matrix or M vector. 460 npts = len(self.x) 461 nq = nr 462 sqrt_alpha = math.sqrt(math.fabs(self.alpha)) 463 if sqrt_alpha<0.0: 464 nq = 0 465 466 err_0 = numpy.zeros([nfunc, nfunc]) 467 c_0 = numpy.zeros(nfunc) 468 nfunc_0 = nfunc 469 nfunc += 1 470 471 a = numpy.zeros([npts+nq, nfunc]) 472 b = numpy.zeros(npts+nq) 473 err = numpy.zeros([nfunc, nfunc]) 474 475 # Construct the a matrix and b vector that represent the problem 476 self._get_matrix(nfunc, nq, a, b) 477 478 c, chi2, rank, n = lstsq(a, b) 479 # Sanity check 480 try: 481 float(chi2) 482 except: 483 chi2 = -1.0 484 self.chi2 = chi2 485 486 inv_cov = numpy.zeros([nfunc,nfunc]) 487 # Get the covariance matrix, defined as inv_cov = a_transposed * a 488 self._get_invcov_matrix(nfunc, nr, a, inv_cov) 489 490 # Compute the reg term size for the output 491 sum_sig, sum_reg = self._get_reg_size(nfunc, nr, a) 492 493 if math.fabs(self.alpha)>0: 494 new_alpha = sum_sig/(sum_reg/self.alpha) 495 else: 496 new_alpha = 0.0 497 self.suggested_alpha = new_alpha 498 499 try: 500 cov = numpy.linalg.pinv(inv_cov) 501 err = math.fabs(chi2/float(npts-nfunc)) * cov 502 except: 503 # We were not able to estimate the errors, 504 # returns an empty covariance matrix 505 print "lstsq:", sys.exc_value 506 print chi2 507 pass 508 509 # Keep a copy of the last output 510 511 print "BACKGROUND =", c[0] 512 self.background = c[0] 513 514 for i in range(nfunc_0): 515 c_0[i] = c[i+1] 516 for j in range(nfunc_0): 517 err_0[i][j] = err[i+1][j+1] 518 519 self.out = c_0 520 self.cov = err_0 521 522 return c_0, err_0 409 523 410 524 def estimate_alpha(self, nfunc): … … 433 547 434 548 # Perform inversion to find the largest alpha 435 out, cov = pr. lstsq(nfunc)549 out, cov = pr.invert(nfunc) 436 550 elapsed = time.time()-starttime 437 551 initial_alpha = pr.alpha … … 440 554 # Try the inversion with the estimated alpha 441 555 pr.alpha = pr.suggested_alpha 442 out, cov = pr. lstsq(nfunc)556 out, cov = pr.invert(nfunc) 443 557 444 558 npeaks = pr.get_peaks(out) … … 458 572 for i in range(10): 459 573 pr.alpha = (0.33)**(i+1)*alpha 460 out, cov = pr. lstsq(nfunc)574 out, cov = pr.invert(nfunc) 461 575 462 576 peaks = pr.get_peaks(out) -
pr_inversion/test/test_output.txt
rb00b487 r9a23253e 5 5 #elapsed=0 6 6 #alpha_estimate=4.58991e-006 7 #C_0=217.876478953+- 0.03949464427868 #C_1=0.30164617344+- 0.4074129227289 #C_2=7.15694574584+- 2.7372648182810 #C_3=-0.79690812244+- 5.718364708311 #C_4=0.948333847403+-1 6.208643264212 #C_5=-0.751612279429+- 28.210889304713 #C_6=-0.0796650783501+-5 2.056830654314 #C_7=-0.518109559543+- 88.606889130915 #C_8=-0.129910831984+-1 41.717830816 #C_9=-0.258134624641+- 215.7674300987 #C_0=217.876478953+-2522.42314999 8 #C_1=0.30164617344+-226.644556603 9 #C_2=7.15694574584+-43.4919208221 10 #C_3=-0.79690812244+-25.8031634616 11 #C_4=0.948333847403+-11.3794625443 12 #C_5=-0.751612279429+-8.64083336005 13 #C_6=-0.0796650783501+-5.73820483268 14 #C_7=-0.518109559543+-3.56660274557 15 #C_8=-0.129910831984+-1.99841316045 16 #C_9=-0.258134624641+-0.756509301101 17 17 <r> <Pr> <dPr> 18 18 0 0 0 19 1.6 22.981 20.513720 3.2 91.986 78.661221 4.8 207.187 164.74122 6.4 368.819 264.12923 8 577.117 359.42724 9.6 832.224 433.16225 11.2 1134.11 470.8926 12.8 1482.46 4 64.81427 14.4 1876.63 419.08628 16 2315.57 360.69329 17.6 2797.81 353.24530 19.2 3321.46 448.35731 20.8 3884.25 611.48132 22.4 4483.6 786.46433 24 5116.72 934.55934 25.6 5780.7 1 030.7835 27.2 6472.66 1 062.2136 28.8 7189.82 1 030.2837 30.4 7929.59 955.68138 32 8689.7 883.9739 33.6 9468.17 878.87140 35.2 10263.4 977.55341 36.8 11073.9 115242 38.4 11898.7 1342.5943 40 12737 1498.344 41.6 13587.7 1587.0845 43.2 14450.2 1596.9546 44.8 15323.5 1537.9947 46.4 16206.5 1445.5848 48 17098 1379.3949 49.6 17996.3 1401.8650 51.2 18899.8 1532.0651 52.8 19806.5 1728.8652 54.4 20714.3 1927.3753 56 21620.8 2073.3954 57.6 22524 2135.2155 59.2 23421.6 2106.9756 60.8 24311.7 2011.2857 62.4 25192.3 1900.7558 64 26061.8 1848.6759 65.6 26918.7 1912.8360 67.2 27761.7 2089.9661 68.8 28589.6 2319.4862 70.4 29401.3 2527.7163 72 30195.6 2658.7664 73.6 30971.4 2685.1765 75.2 31727.2 2611.6866 76.8 32461.5 2478.0567 78.4 33172.7 2357.1968 80 33858.8 2333.8169 81.6 34517.8 2453.4370 83.2 35147.6 2684.6171 84.8 35745.9 2945.1372 86.4 36310.5 3152.1873 88 36839.5 3249.5974 89.6 37330.9 3217.0475 91.2 37783.1 3074.6176 92.8 38194.6 2886.0977 94.4 38564.3 2752.5878 96 38891.1 2773.0379 97.6 39174.1 2973.0280 99.2 39412.6 3281.6281 100.8 39605.8 3587.5982 102.4 39752.6 3795.9583 104 39851.9 3850.5684 105.6 39902.4 3741.3285 107.2 39902.5 3510.0386 108.8 39850.1 3255.5887 110.4 39743.2 3120.2588 112 39579.6 3218.6189 113.6 39356.7 3539.2690 115.2 39072.3 3954.9291 116.8 38724.1 4317.7392 118.4 38310.3 4517.0893 120 37829 4494.8694 121.6 37278.9 4251.4495 123.2 36659 3856.5896 124.8 35968.4 3465.7897 126.4 35206.6 3306.2298 128 34373 3535.9799 129.6 33466.9 4074.35100 131.2 32487.5 4693.64101 132.8 31433.4 5186.14102 134.4 30303.1 5411.6103 136 29094.1 5295.75104 137.6 27803.8 4830.98105 139.2 26428.9 4092.02106 140.8 24965.9 3287.78107 142.4 23411.1 2858.06108 144 21761.3 3246.41109 145.6 20013.4 4237.67110 147.2 18165.4 5345.59111 148.8 16216.4 6256.27112 150.4 14166.9 6786.26113 152 12019 6829.04114 153.6 9776.61 6338.91115 155.2 7445.34 5326.31116 156.8 5032.57 3854117 158.4 2547.22 2030.419 1.6 22.981 7.48836 20 3.2 91.986 29.737 21 4.8 207.187 66.1302 22 6.4 368.819 115.749 23 8 577.117 177.508 24 9.6 832.224 250.315 25 11.2 1134.11 333.205 26 12.8 1482.46 425.442 27 14.4 1876.63 526.549 28 16 2315.57 636.272 29 17.6 2797.81 754.499 30 19.2 3321.46 881.145 31 20.8 3884.25 1016.06 32 22.4 4483.6 1158.97 33 24 5116.72 1309.47 34 25.6 5780.7 1467.05 35 27.2 6472.66 1631.18 36 28.8 7189.82 1801.33 37 30.4 7929.59 1977.02 38 32 8689.7 2157.88 39 33.6 9468.17 2343.55 40 35.2 10263.4 2533.72 41 36.8 11073.9 2728.06 42 38.4 11898.7 2926.2 43 40 12737 3127.75 44 41.6 13587.7 3332.27 45 43.2 14450.2 3539.34 46 44.8 15323.5 3748.54 47 46.4 16206.5 3959.49 48 48 17098 4171.83 49 49.6 17996.3 4385.18 50 51.2 18899.8 4599.19 51 52.8 19806.5 4813.46 52 54.4 20714.3 5027.56 53 56 21620.8 5241.07 54 57.6 22524 5453.53 55 59.2 23421.6 5664.53 56 60.8 24311.7 5873.68 57 62.4 25192.3 6080.65 58 64 26061.8 6285.13 59 65.6 26918.7 6486.85 60 67.2 27761.7 6685.53 61 68.8 28589.6 6880.91 62 70.4 29401.3 7072.73 63 72 30195.6 7260.72 64 73.6 30971.4 7444.63 65 75.2 31727.2 7624.23 66 76.8 32461.5 7799.3 67 78.4 33172.7 7969.61 68 80 33858.8 8134.91 69 81.6 34517.8 8294.9 70 83.2 35147.6 8449.24 71 84.8 35745.9 8597.54 72 86.4 36310.5 8739.34 73 88 36839.5 8874.2 74 89.6 37330.9 9001.64 75 91.2 37783.1 9121.2 76 92.8 38194.6 9232.39 77 94.4 38564.3 9334.72 78 96 38891.1 9427.68 79 97.6 39174.1 9510.74 80 99.2 39412.6 9583.35 81 100.8 39605.8 9644.98 82 102.4 39752.6 9695.14 83 104 39851.9 9733.38 84 105.6 39902.4 9759.36 85 107.2 39902.5 9772.74 86 108.8 39850.1 9773.24 87 110.4 39743.2 9760.52 88 112 39579.6 9734.22 89 113.6 39356.7 9693.89 90 115.2 39072.3 9639.07 91 116.8 38724.1 9569.28 92 118.4 38310.3 9484.1 93 120 37829 9383.18 94 121.6 37278.9 9266.24 95 123.2 36659 9133 96 124.8 35968.4 8983.13 97 126.4 35206.6 8816.15 98 128 34373 8631.45 99 129.6 33466.9 8428.29 100 131.2 32487.5 8205.96 101 132.8 31433.4 7963.92 102 134.4 30303.1 7701.95 103 136 29094.1 7420.23 104 137.6 27803.8 7119.29 105 139.2 26428.9 6799.69 106 140.8 24965.9 6461.63 107 142.4 23411.1 6104.48 108 144 21761.3 5726.35 109 145.6 20013.4 5323.88 110 147.2 18165.4 4892.46 111 148.8 16216.4 4426.74 112 150.4 14166.9 3921.46 113 152 12019 3372.52 114 153.6 9776.61 2777.82 115 155.2 7445.34 2138.05 116 156.8 5032.57 1456.94 117 158.4 2547.22 741.175 -
pr_inversion/test/utest_invertor.py
rb00b487 r9a23253e 63 63 self.x_in[i] = 1.0*(i+1) 64 64 65 def test_has_bck_flag(self): 66 """ 67 Tests the has_bck flag operations 68 """ 69 self.assertEqual(self.invertor.has_bck, False) 70 self.invertor.has_bck=True 71 self.assertEqual(self.invertor.has_bck, True) 72 def doit_float(): 73 self.invertor.has_bck = 2.0 74 def doit_str(): 75 self.invertor.has_bck = 'a' 76 77 self.assertRaises(ValueError, doit_float) 78 self.assertRaises(ValueError, doit_str) 79 65 80 66 81 def testset_dmax(self): … … 177 192 # Perform inversion 178 193 out, cov = self.invertor.invert_optimize(10) 194 #out, cov = self.invertor.invert(10) 179 195 # This is a very specific case 180 196 # We should make sure it always passes
Note: See TracChangeset
for help on using the changeset viewer.