Changeset 7708 for trunk/projects/slimcurve/src/main/c/EcfSingle.c
 Timestamp:
 05/20/11 22:24:58 (9 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/projects/slimcurve/src/main/c/EcfSingle.c
r7706 r7708 4 4 5 5 This file contains functions for single transient analysis. 6 Utility code is found in EcfUtil.c. 6 Utility code is found in EcfUtil.c and global analysis code in 7 EcfGlobal.c. 7 8 */ 8 9 … … 13 14 14 15 /* Predeclarations */ 16 int GCI_marquardt_step(float x[], float y[], int ndata, 17 noise_type noise, float sig[], 18 float param[], int paramfree[], int nparam, 19 restrain_type restrain, 20 void (*fitfunc)(float, float [], float *, float [], int), 21 float yfit[], float dy[], 22 float **covar, float **alpha, float *chisq, 23 float *alambda, int *pmfit, float *pochisq, float *paramtry, float *beta, float *dparam); 15 24 int GCI_marquardt_step_instr(float xincr, float y[], 16 25 int ndata, int fit_start, int fit_end, … … 22 31 float yfit[], float dy[], 23 32 float **covar, float **alpha, float *chisq, 24 float *alambda, int *pmfit, float *pochisq, float *paramtry, float *beta, float *dparam, 33 float *alambda, int *pmfit, float *pochisq, float *paramtry, float *beta, float *dparam, 25 34 float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv, 26 35 int *pfnvals_len, int *pdy_dparam_nparam_size); 36 int GCI_marquardt_estimate_errors(float **alpha, int nparam, int mfit, 37 float d[], float **v, float interval); 27 38 28 39 /* Globals */ … … 349 360 float *chisq, float chisq_target) 350 361 { 351 int tries=1, division=3; // the data 362 int tries=1, division=3; // the data 352 363 float local_chisq=3.0e38, oldChisq=3.0e38, oldZ, oldA, oldTau, *validFittedArray; // local_chisq a very high float but below oldChisq 353 364 354 365 if (fitted==NULL) // we require chisq but have not supplied a "fitted" array so must malloc one 355 366 { … … 357 368 } 358 369 else validFittedArray = fitted; 359 370 360 371 if (instr==NULL) // no instrument/prompt has been supplied 361 372 { … … 363 374 Z, A, tau, validFittedArray, residuals, &local_chisq, division); 364 375 365 while (local_chisq>chisq_target && (local_chisq<=oldChisq) && tries<MAXREFITS) 376 while (local_chisq>chisq_target && (local_chisq<=oldChisq) && tries<MAXREFITS) 366 377 { 367 378 oldChisq = local_chisq; … … 369 380 oldA = *A; 370 381 oldTau = *tau; 371 // division++; 372 division+=division/3; 382 // division++; 383 division+=division/3; 373 384 tries++; 374 385 GCI_triple_integral(xincr, y, fit_start, fit_end, noise, sig, 375 386 Z, A, tau, validFittedArray, residuals, &local_chisq, division); 376 } 387 } 377 388 } 378 389 else … … 381 392 Z, A, tau, validFittedArray, residuals, &local_chisq, division); 382 393 383 while (local_chisq>chisq_target && (local_chisq<=oldChisq) && tries<MAXREFITS) 394 while (local_chisq>chisq_target && (local_chisq<=oldChisq) && tries<MAXREFITS) 384 395 { 385 oldChisq = local_chisq; 396 oldChisq = local_chisq; 386 397 oldZ = *Z; 387 398 oldA = *A; 388 399 oldTau = *tau; 389 // division++; 400 // division++; 390 401 division+=division/3; 391 402 tries++; … … 393 404 Z, A, tau, validFittedArray, residuals, &local_chisq, division); 394 405 395 } 406 } 396 407 } 397 408 … … 403 414 *tau = oldTau; 404 415 } 405 416 406 417 if (chisq!=NULL) *chisq = local_chisq; 407 418 408 if (fitted==NULL) 419 if (fitted==NULL) 409 420 { 410 421 free (validFittedArray); … … 454 465 noise=NOISE_POISSON_FIT, the variances are taken to be 455 466 max(yfit[i],15). If noise=NOISE_GAUSSIAN_FIT, the variances are taken to be 456 yfit[i] and if noise=NOISE_MLE then a optimisation is for the maximum 457 likelihood approximation (Laurence and Chromy in press 2010 and 467 yfit[i] and if noise=NOISE_MLE then a optimisation is for the maximum 468 likelihood approximation (Laurence and Chromy in press 2010 and 458 469 G. Nishimura, and M. Tamura Phys Med Biol 2005). 459 470 460 471 The input array paramfree[0..nparam1] indicates 461 472 by nonzero entries those components that should be held fixed at … … 488 499 */ 489 500 490 /* This functions does the whole job */ 491 501 /* These two functions do the whole job */ 492 502 int GCI_marquardt(float x[], float y[], int ndata, 493 503 noise_type noise, float sig[], … … 499 509 float chisq_delta, float chisq_percent, float **erraxes) 500 510 { 501 // Empty fn to allow compile in TRI2 511 float alambda, ochisq; 512 int mfit; 513 float evals[MAXFIT]; 514 int i, k, itst, itst_max; 515 516 float ochisq2, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT]; 517 518 itst_max = (restrain == ECF_RESTRAIN_DEFAULT) ? 4 : 6; 519 520 mfit = 0; 521 for (i=0; i<nparam; i++) { 522 if (paramfree[i]) mfit++; 523 } 524 525 alambda = 1; 526 if (GCI_marquardt_step(x, y, ndata, noise, sig, 527 param, paramfree, nparam, restrain, 528 fitfunc, fitted, residuals, 529 covar, alpha, chisq, &alambda, 530 &mfit, &ochisq2, paramtry, beta, dparam) != 0) { 531 return 1; 532 } 533 534 k = 1; /* Iteration counter */ 535 itst = 0; 536 for (;;) { 537 k++; 538 if (k > MAXITERS) { 539 return 2; 540 } 541 542 ochisq = *chisq; 543 if (GCI_marquardt_step(x, y, ndata, noise, sig, 544 param, paramfree, nparam, restrain, 545 fitfunc, fitted, residuals, 546 covar, alpha, chisq, &alambda, 547 &mfit, &ochisq2, paramtry, beta, dparam) != 0) { 548 return 3; 549 } 550 551 if (*chisq > ochisq) 552 itst = 0; 553 else if (ochisq  *chisq < chisq_delta) 554 itst++; 555 556 if (itst < itst_max) continue; 557 558 /* Endgame */ 559 alambda = 0.0; 560 if (GCI_marquardt_step(x, y, ndata, noise, sig, 561 param, paramfree, nparam, restrain, 562 fitfunc, fitted, residuals, 563 covar, alpha, chisq, &alambda, 564 &mfit, &ochisq2, paramtry, beta, dparam) != 0) { 565 return 4; 566 } 567 568 if (erraxes == NULL) 569 return k; 570 571 //TODO ARG 572 //if (GCI_marquardt_estimate_errors(alpha, nparam, mfit, evals, 573 // erraxes, chisq_percent) != 0) { 574 // return 5; 575 //} 576 577 break; /* We're done now */ 578 } 579 580 return k; 502 581 } 503 582 … … 529 608 float *fnvals=NULL, **dy_dparam_pure=NULL, **dy_dparam_conv=NULL; 530 609 int fnvals_len=0, dy_dparam_nparam_size=0; 531 float ochisq2, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT]; 610 float ochisq2, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT]; 532 611 533 612 itst_max = (restrain == ECF_RESTRAIN_DEFAULT) ? 4 : 6; … … 545 624 param, paramfree, nparam, restrain, 546 625 fitfunc, fitted, residuals, 547 covar, alpha, chisq, &alambda, 626 covar, alpha, chisq, &alambda, 548 627 &mfit2, &ochisq2, paramtry, beta, dparam, 549 628 &fnvals, &dy_dparam_pure, &dy_dparam_conv, … … 569 648 param, paramfree, nparam, restrain, 570 649 fitfunc, fitted, residuals, 571 covar, alpha, chisq, &alambda, 650 covar, alpha, chisq, &alambda, 572 651 &mfit2, &ochisq2, paramtry, beta, dparam, 573 652 &fnvals, &dy_dparam_pure, &dy_dparam_conv, … … 592 671 param, paramfree, nparam, restrain, 593 672 fitfunc, fitted, residuals, 594 covar, alpha, chisq, &alambda, 673 covar, alpha, chisq, &alambda, 595 674 &mfit2, &ochisq2, paramtry, beta, dparam, 596 675 &fnvals, &dy_dparam_pure, &dy_dparam_conv, … … 605 684 } 606 685 686 //TODO ARG this estimate errors call was deleted in my latest version 687 // if (GCI_marquardt_estimate_errors(alpha, nparam, mfit, evals, 688 // erraxes, chisq_percent) != 0) { 689 // do_frees 690 // return 5; 691 // } 692 607 693 break; /* We're done now */ 608 694 } … … 612 698 } 613 699 #undef do_frees 700 701 702 //TODO ARG deleted in my latest version 703 int GCI_marquardt_step(float x[], float y[], int ndata, 704 noise_type noise, float sig[], 705 float param[], int paramfree[], int nparam, 706 restrain_type restrain, 707 void (*fitfunc)(float, float [], float *, float [], int), 708 float yfit[], float dy[], 709 float **covar, float **alpha, float *chisq, 710 float *alambda, int *pmfit, float *pochisq, float *paramtry, float *beta, float *dparam) 711 { 712 int j, k, l, ret; 713 // static int mfit; // was static but now thread safe 714 // static float ochisq, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT]; // was static but now thread safe 715 int mfit = *pmfit; 716 float ochisq = *pochisq; 717 718 if (nparam > MAXFIT) 719 return 1; 720 721 /* Initialisation */ 722 /* We assume we're given sensible starting values for param[] */ 723 if (*alambda < 0.0) { 724 for (mfit=0, j=0; j<nparam; j++) 725 if (paramfree[j]) 726 mfit++; 727 728 if (GCI_marquardt_compute_fn(x, y, ndata, noise, sig, 729 param, paramfree, nparam, fitfunc, 730 yfit, dy, 731 alpha, beta, chisq, 0.0, *alambda) != 0) 732 return 2; 733 734 *alambda = 0.001; 735 ochisq = (*chisq); 736 for (j=0; j<nparam; j++) 737 paramtry[j] = param[j]; 738 } 739 740 /* Alter linearised fitting matrix by augmenting diagonal elements */ 741 for (j=0; j<mfit; j++) { 742 for (k=0; k<mfit; k++) 743 covar[j][k] = alpha[j][k]; 744 745 covar[j][j] = alpha[j][j] * (1.0 + (*alambda)); 746 dparam[j] = beta[j]; 747 } 748 749 /* Matrix solution; GCI_solve solves Ax=b rather than AX=B */ 750 if (GCI_solve(covar, mfit, dparam) != 0) 751 return 1; 752 753 //TODO ARG GCI_gauss_jordan would invert the covar matrix as a side effect 754 /* Once converged, evaluate covariance matrix */ 755 if (*alambda == 0) { 756 if (GCI_marquardt_compute_fn_final(x, y, ndata, noise, sig, 757 param, paramfree, nparam, fitfunc, 758 yfit, dy, chisq) != 0) 759 return 3; 760 if (mfit < nparam) { /* no need to do this otherwise */ 761 GCI_covar_sort(covar, nparam, paramfree, mfit); 762 GCI_covar_sort(alpha, nparam, paramfree, mfit); 763 } 764 return 0; 765 } 766 767 /* Did the trial succeed? */ 768 for (j=0, l=0; l<nparam; l++) 769 if (paramfree[l]) 770 paramtry[l] = param[l] + dparam[j++]; 771 772 if (restrain == ECF_RESTRAIN_DEFAULT) 773 ret = check_ecf_params (paramtry, nparam, fitfunc); 774 else 775 ret = check_ecf_user_params (paramtry, nparam, fitfunc); 776 777 if (ret != 0) { 778 /* Bad parameters, increase alambda and return */ 779 *alambda *= 10.0; 780 return 0; 781 } 782 783 if (GCI_marquardt_compute_fn(x, y, ndata, noise, sig, 784 paramtry, paramfree, nparam, fitfunc, 785 yfit, dy, covar, dparam, 786 chisq, ochisq, *alambda) != 0) 787 return 2; 788 789 /* Success, accept the new solution */ 790 if (*chisq < ochisq) { 791 *alambda *= 0.1; 792 ochisq = *chisq; 793 for (j=0; j<mfit; j++) { 794 for (k=0; k<mfit; k++) 795 alpha[j][k] = covar[j][k]; 796 beta[j] = dparam[j]; 797 } 798 for (l=0; l<nparam; l++) 799 param[l] = paramtry[l]; 800 } else { /* Failure, increase alambda and return */ 801 *alambda *= 10.0; 802 *chisq = ochisq; 803 } 804 805 return 0; 806 } 614 807 615 808 … … 623 816 float yfit[], float dy[], 624 817 float **covar, float **alpha, float *chisq, 625 float *alambda, int *pmfit, float *pochisq, float *paramtry, float *beta, float *dparam, 818 float *alambda, int *pmfit, float *pochisq, float *paramtry, float *beta, float *dparam, 626 819 float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv, 627 820 int *pfnvals_len, int *pdy_dparam_nparam_size) … … 631 824 // static float ochisq, paramtry[MAXFIT], beta[MAXFIT], dparam[MAXFIT]; // was static but now thread safe 632 825 633 /*if (nparam > MAXFIT) { 634 printf("GCI_marq_step_instr returns 10\n"); 826 if (nparam > MAXFIT) 635 827 return 10; 636 } 637 if (xincr <= 0) { 638 printf("GCI_marq_step_instr returns 11\n"); 828 if (xincr <= 0) 639 829 return 11; 640 } 641 if (fit_start < 0  fit_start > fit_end  fit_end > ndata) { 642 printf("GCI_marq_step_instr returns 12\n"); 830 if (fit_start < 0  fit_start > fit_end  fit_end > ndata) 643 831 return 12; 644 }*/645 832 646 833 /* Initialisation */ … … 656 843 param, paramfree, nparam, fitfunc, 657 844 yfit, dy, alpha, beta, chisq, 0.0, 658 *alambda, 845 *alambda, 659 846 pfnvals, pdy_dparam_pure, pdy_dparam_conv, 660 pfnvals_len, pdy_dparam_nparam_size) != 0) {847 pfnvals_len, pdy_dparam_nparam_size) != 0) 661 848 return 2; 662 }663 849 664 850 *alambda = 0.001; … … 678 864 } 679 865 680 /* Matrix solution; GCI_ solvesolves Ax=b rather than AX=B */681 if (GCI_solve(covar, *pmfit, dparam) != 0) {866 /* Matrix solution; GCI_gauss_jordan solves Ax=b rather than AX=B */ 867 if (GCI_solve(covar, *pmfit, dparam) != 0) 682 868 return 1; 683 } 684 685 //TODO need to make sure covar gets inverted. Previously the Gauss 686 // Jordan solution would invert covar as a side effect. 687 869 870 //TODO ARG covar needs to get inverted; previously inverted as a side effect 688 871 /* Once converged, evaluate covariance matrix */ 689 872 if (*alambda == 0) { … … 692 875 instr, ninstr, noise, sig, 693 876 param, paramfree, nparam, fitfunc, 694 yfit, dy, chisq, 877 yfit, dy, chisq, 695 878 pfnvals, pdy_dparam_pure, pdy_dparam_conv, 696 pfnvals_len, pdy_dparam_nparam_size) != 0) {879 pfnvals_len, pdy_dparam_nparam_size) != 0) 697 880 return 3; 698 }699 881 if (*pmfit < nparam) { /* no need to do this otherwise */ 700 882 GCI_covar_sort(covar, nparam, paramfree, *pmfit); … … 704 886 } 705 887 706 /* Did the trial succeed? * 707 int allzeroes = 1; 888 //TODO ARG c/b special case if dparam all zeroes; don't have to recalc alpha & beta 889 /* Did the trial succeed? */ 708 890 for (j=0, l=0; l<nparam; l++) 709 if (paramfree[l]) { 710 //TODO I don't think it ever recovers if all dparams are zero; alambda grows and grows, goes from int to inf!! 711 // could be an early way to detect we're done. 712 // at least, if delta params is all zeroes, there is no point in calculating a new alpha and beta. 713 // printf("paramtry[%d] was %f becomes %f, param[%d] or %f plus dparam[%d] or %f\n", l, paramtry[l], param[l] + dparam[j], l, param[l], j, dparam[j]); 714 if (dparam[j] != 0.0) allzeroes = 0; 891 if (paramfree[l]) 715 892 paramtry[l] = param[l] + dparam[j++]; 716 }717 //TODO experimental:718 // gave about a 10% speedup719 if (allzeroes) return 12345; */720 893 721 894 if (restrain == ECF_RESTRAIN_DEFAULT) … … 729 902 return 0; 730 903 } 904 731 905 if (GCI_marquardt_compute_fn_instr(xincr, y, ndata, fit_start, fit_end, 732 906 instr, ninstr, noise, sig, … … 735 909 chisq, *pochisq, *alambda, 736 910 pfnvals, pdy_dparam_pure, pdy_dparam_conv, 737 pfnvals_len, pdy_dparam_nparam_size) != 0) {911 pfnvals_len, pdy_dparam_nparam_size) != 0) 738 912 return 2; 739 }740 913 741 914 /* Success, accept the new solution */ 742 915 if (*chisq < *pochisq) { 743 // printf("success, alambda goes from %f to %f\n", *alambda, (*alambda)*0.1);744 916 *alambda *= 0.1; 745 917 *pochisq = *chisq; … … 752 924 param[l] = paramtry[l]; 753 925 } else { /* Failure, increase alambda and return */ 754 //TODO if alambda goes to inf, can we recover?755 // printf("failure, alambda goes from %f to %f\n", *alambda, (*alambda)*10.0);756 926 *alambda *= 10.0; 757 927 *chisq = *pochisq; … … 789 959 determined the fitted convolved response. 790 960 791 This is the variant which handles an instrument response. */ 792 961 Again there are two variants of this function, corresponding to the 962 two variants of the Marquardt function. 963 */ 964 965 //TODO ARG deleted in my version; needs to get deNRed!!!! 966 int GCI_marquardt_compute_fn(float x[], float y[], int ndata, 967 noise_type noise, float sig[], 968 float param[], int paramfree[], int nparam, 969 void (*fitfunc)(float, float [], float *, float [], int), 970 float yfit[], float dy[], 971 float **alpha, float beta[], float *chisq, float old_chisq, 972 float alambda) 973 { 974 int i, j, k, l, m, mfit; 975 float wt, sig2i, y_ymod, dy_dparam[MAXFIT]; 976 977 for (j=0, mfit=0; j<nparam; j++) 978 if (paramfree[j]) 979 mfit++; 980 981 //TODO ARG having this section { } of code here is just temporary; 982 // o'wise have to define these vars at top of function 983 { 984 float alpha_weight[MAXBINS]; 985 float beta_weight[MAXBINS]; 986 int q; 987 float weight; 988 989 int i_free; 990 int j_free; 991 float dot_product; 992 float beta_sum; 993 float dy_dparam_k_i; 994 995 *chisq = 0.0f; 996 997 switch (noise) { 998 case NOISE_CONST: 999 { 1000 for (q = 0; q < ndata; ++q) { 1001 (*fitfunc)(x[q], param, &yfit[q], dy_dparam, nparam); 1002 yfit[q] += param[0]; 1003 dy[q] = y[q]  yfit[q]; 1004 weight = 1.0f; //TODO ARG this should be 1.0f / sig[0] ! 1005 alpha_weight[q] = weight; // 1 1006 weight *= dy[q]; 1007 beta_weight[q] = weight; // dy[q] 1008 weight *= dy[q]; 1009 *chisq += weight; // dy[q] * dy[q] 1010 } 1011 break; 1012 } 1013 case NOISE_GIVEN: 1014 { 1015 for (q = 0; q < ndata; ++q) { 1016 (*fitfunc)(x[q], param, &yfit[q], dy_dparam, nparam); 1017 yfit[q] += param[0]; 1018 dy[q] = y[q]  yfit[q]; 1019 weight = 1.0f / (sig[q] * sig[q]); 1020 alpha_weight[q] = weight; // 1 / (sig[q] * sig[q]) 1021 weight *= dy[q]; 1022 beta_weight[q] = weight; // dy[q] / (sig[q] * sig[q]) 1023 weight *= dy[q]; 1024 *chisq += weight; // (dy[q] * dy[q]) / (sig[q] * sig[q]) 1025 } 1026 break; 1027 } 1028 case NOISE_POISSON_DATA: 1029 { 1030 for (q = 0; q < ndata; ++q) { 1031 (*fitfunc)(x[q], param, &yfit[q], dy_dparam, nparam); 1032 yfit[q] += param[0]; 1033 dy[q] = y[q]  yfit[q]; 1034 weight = (y[q] > 15 ? 1.0f / y[q] : 1.0f / 15); 1035 alpha_weight[q] = weight; // 1 / sig(q) 1036 weight *= dy[q]; 1037 beta_weight[q] = weight; // dy[q] / sig(q) 1038 weight *= dy[q]; 1039 *chisq += weight; // (dy[q] * dy[q]) / sig(q) 1040 } 1041 break; 1042 } 1043 case NOISE_POISSON_FIT: 1044 { 1045 for (q = 0; q < ndata; ++q) { 1046 (*fitfunc)(x[q], param, &yfit[q], dy_dparam, nparam); 1047 yfit[q] += param[0]; 1048 dy[q] = y[q]  yfit[q]; 1049 weight = (yfit[q] > 15 ? 1.0f / yfit[q] : 1.0f / 15); 1050 alpha_weight[q] = weight; // 1 / sig(q) 1051 weight *= dy[q]; 1052 beta_weight[q] = weight; // dy(q) / sig(q) 1053 weight *= dy[q]; 1054 *chisq += weight; // (dy(q) * dy(q)) / sig(q) 1055 } 1056 break; 1057 } 1058 case NOISE_GAUSSIAN_FIT: 1059 { 1060 for (q = 0; q < ndata; ++q) { 1061 (*fitfunc)(x[q], param, &yfit[q], dy_dparam, nparam); 1062 yfit[q] += param[0]; 1063 dy[q] = y[q]  yfit[q]; 1064 weight = (yfit[q] > 1.0f ? 1.0f / yfit[q] : 1.0f); 1065 alpha_weight[q] = weight; // 1 / sig(q) 1066 weight *= dy[q]; 1067 beta_weight[q] = weight; // dy[q] / sig(q) 1068 weight *= dy[q]; 1069 *chisq += weight; 1070 } 1071 break; 1072 } 1073 case NOISE_MLE: 1074 { 1075 for (q = 0; q < ndata; ++q) { 1076 (*fitfunc)(x[q], param, &yfit[q], dy_dparam, nparam); 1077 yfit[q] += param[0]; 1078 dy[q] = y[q]  yfit[q]; 1079 weight = (yfit[q] > 1 ? 1.0f / yfit[i] : 1.0f); 1080 alpha_weight[q] = weight * y[q] / yfit[q]; //TODO was y_ymod = y[i]/yfit[i], but y_ymod was float, s/b modulus? 1081 beta_weight[q] = dy[q] * weight; 1082 *chisq += (0.0f == y[q]) 1083 ? 2.0 * yfit[q] 1084 : 2.0 * (yfit[q]  y[q])  2.0 * y[q] * log(yfit[q] / y[q]); 1085 } 1086 if (*chisq <= 0.0f) { 1087 *chisq = 1.0e38f; // don't let chisq=0 through yfit being all ve 1088 } 1089 break; 1090 } 1091 default: 1092 { 1093 return 3; 1094 } 1095 } 1096 1097 // Check if chi square worsened: 1098 if (0.0f != old_chisq && *chisq >= old_chisq) { 1099 // don't bother to set up the matrices for solution 1100 return 0; 1101 } 1102 1103 i_free = 0; 1104 // for all columns 1105 for (i = 0; i < nparam; ++i) { 1106 if (paramfree[i]) { 1107 j_free = 0; 1108 beta_sum = 0.0f; 1109 // row loop, only need to consider lower triangle 1110 for (j = 0; j <= i; ++j) { 1111 if (paramfree[j]) { 1112 dot_product = 0.0f; 1113 if (0 == j_free) { 1114 // for all data 1115 for (k = 0; k < ndata; ++k) { 1116 //TODO ARG just to get this to compile, for now: 1117 //TODO ARG from the _instr version: dy_dparam_k_i = (*pdy_dparam_conv)[k][i]; 1118 //TODO ARG from the instr version dot_product += dy_dparam_k_i * (*pdy_dparam_conv)[k][j] * alpha_weight[k]; //TODO make it [i][k] and just *dy_dparam++ it. 1119 //TODO ARG from the instr version beta_sum += dy_dparam_k_i * beta_weight[k]; 1120 } 1121 } 1122 else { 1123 // for all data 1124 for (k = 0; k < ndata; ++k) { 1125 //TODO ARG from the instr version: dot_product += (*pdy_dparam_conv)[k][i] * (*pdy_dparam_conv)[k][j] * alpha_weight[k]; 1126 } 1127 } // k loop 1128 1129 alpha[j_free][i_free] = alpha[i_free][j_free] = dot_product; 1130 // if (i_free != j_free) { 1131 // // matrix is symmetric 1132 // alpha[i_free][j_free] = dot_product; //TODO dotProduct s/b including fixed parameters????!!! 1133 // } 1134 ++j_free; 1135 } 1136 } // j loop 1137 beta[i_free] = beta_sum; 1138 ++i_free; 1139 } 1140 //else printf("param not free %d\n", i); 1141 } // i loop 1142 } 1143 1144 return 0; 1145 1146 1147 /* Initialise (symmetric) alpha, beta */ 1148 //TODO ARG FRI 1149 //for (j=0; j<mfit; j++) { 1150 // for (k=0; k<=j; k++) 1151 // alpha[j][k] = 0.0; 1152 // beta[j] = 0.0; 1153 //} 1154 1155 /* Calculation of the fitting data will depend upon the type of 1156 noise and the type of instrument response */ 1157 1158 /* There's no convolution involved in this function. This is then 1159 fairly straightforward, depending only upon the type of noise 1160 present. Since we calculate the standard deviation at every 1161 point in a different way depending upon the type of noise, we 1162 will place our switch(noise) around the entire loop. */ 1163 1164 switch (noise) { 1165 case NOISE_CONST: 1166 *chisq = 0.0; 1167 /* Summation loop over all data */ 1168 for (i=0; i<ndata; i++) { 1169 (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 1170 yfit[i] += param[0]; 1171 dy_dparam[0] = 1; 1172 dy[i] = y[i]  yfit[i]; 1173 for (j=0, l=0; l<nparam; l++) { 1174 if (paramfree[l]) { 1175 wt = dy_dparam[l]; /* taken away the *sig2i from here */ 1176 for (k=0, m=0; m<=l; m++) 1177 if (paramfree[m]) 1178 alpha[j][k++] += wt * dy_dparam[m]; 1179 beta[j] += dy[i] * wt; 1180 j++; 1181 } 1182 } 1183 /* And find chi^2 */ 1184 *chisq += dy[i] * dy[i]; 1185 } 1186 1187 /* Now divide everything by sigma^2 */ 1188 sig2i = 1.0 / (sig[0] * sig[0]); 1189 *chisq *= sig2i; 1190 for (j=0; j<mfit; j++) { 1191 for (k=0; k<=j; k++) 1192 alpha[j][k] *= sig2i; 1193 beta[j] *= sig2i; 1194 } 1195 break; 1196 1197 case NOISE_GIVEN: /* This is essentially the NR version */ 1198 *chisq = 0.0; 1199 /* Summation loop over all data */ 1200 for (i=0; i<ndata; i++) { 1201 (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 1202 yfit[i] += param[0]; 1203 dy_dparam[0] = 1; 1204 sig2i = 1.0 / (sig[i] * sig[i]); 1205 dy[i] = y[i]  yfit[i]; 1206 for (j=0, l=0; l<nparam; l++) { 1207 if (paramfree[l]) { 1208 wt = dy_dparam[l] * sig2i; 1209 for (k=0, m=0; m<=l; m++) 1210 if (paramfree[m]) 1211 alpha[j][k++] += wt * dy_dparam[m]; 1212 beta[j] += wt * dy[i]; 1213 j++; 1214 } 1215 } 1216 /* And find chi^2 */ 1217 *chisq += dy[i] * dy[i] * sig2i; 1218 } 1219 break; 1220 1221 case NOISE_POISSON_DATA: 1222 *chisq = 0.0; 1223 /* Summation loop over all data */ 1224 for (i=0; i<ndata; i++) { 1225 (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 1226 yfit[i] += param[0]; 1227 dy_dparam[0] = 1; 1228 sig2i = (y[i] > 15 ? 1.0/y[i] : 1.0/15); 1229 dy[i] = y[i]  yfit[i]; 1230 for (j=0, l=0; l<nparam; l++) { 1231 if (paramfree[l]) { 1232 wt = dy_dparam[l] * sig2i; 1233 for (k=0, m=0; m<=l; m++) 1234 if (paramfree[m]) 1235 alpha[j][k++] += wt * dy_dparam[m]; 1236 beta[j] += wt * dy[i]; 1237 j++; 1238 } 1239 } 1240 /* And find chi^2 */ 1241 *chisq += dy[i] * dy[i] * sig2i; 1242 } 1243 break; 1244 1245 case NOISE_POISSON_FIT: 1246 *chisq = 0.0; 1247 // Summation loop over all data 1248 for (i=0; i<ndata; i++) { 1249 (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 1250 yfit[i] += param[0]; 1251 dy_dparam[0] = 1; 1252 sig2i = (yfit[i] > 15 ? 1.0/yfit[i] : 1.0/15); 1253 dy[i] = y[i]  yfit[i]; 1254 for (j=0, l=0; l<nparam; l++) { 1255 if (paramfree[l]) { 1256 wt = dy_dparam[l] * sig2i; 1257 for (k=0, m=0; m<=l; m++) 1258 if (paramfree[m]) 1259 alpha[j][k++] += wt * dy_dparam[m]; 1260 beta[j] += wt * dy[i]; 1261 j++; 1262 } 1263 } 1264 // And find chi^2 1265 *chisq += dy[i] * dy[i] * sig2i; 1266 } 1267 break; 1268 1269 case NOISE_GAUSSIAN_FIT: 1270 *chisq = 0.0; 1271 // Summation loop over all data 1272 for (i=0; i<ndata; i++) { 1273 (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 1274 yfit[i] += param[0]; 1275 dy_dparam[0] = 1; 1276 sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); 1277 dy[i] = y[i]  yfit[i]; 1278 for (j=0, l=0; l<nparam; l++) { 1279 if (paramfree[l]) { 1280 wt = dy_dparam[l] * sig2i; 1281 for (k=0, m=0; m<=l; m++) 1282 if (paramfree[m]) 1283 alpha[j][k++] += wt * dy_dparam[m]; 1284 beta[j] += wt * dy[i]; 1285 j++; 1286 } 1287 } 1288 // And find chi^2 1289 *chisq += dy[i] * dy[i] * sig2i; 1290 } 1291 break; 1292 1293 case NOISE_MLE: 1294 *chisq = 0.0; 1295 /* Summation loop over all data */ 1296 for (i=0; i<ndata; i++) { 1297 (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 1298 yfit[i] += param[0]; 1299 dy_dparam[0] = 1; 1300 sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); 1301 dy[i] = y[i]  yfit[i]; 1302 y_ymod=y[i]/yfit[i]; 1303 for (j=0, l=0; l<nparam; l++) { 1304 if (paramfree[l]) { 1305 wt = dy_dparam[l] * sig2i; 1306 for (k=0, m=0; m<=l; m++) 1307 if (paramfree[m]) 1308 alpha[j][k++] += wt*dy_dparam[m]*y_ymod; //wt * dy_dparam[m]; 1309 beta[j] += wt * dy[i]; 1310 j++; 1311 } 1312 } 1313 // And find chi^2 1314 if (yfit[i]<=0.0) 1315 ; // do nothing 1316 else if (y[i]==0.0) 1317 *chisq += 2.0*yfit[i]; // to avoid NaN from log 1318 else 1319 *chisq += 2.0*(yfit[i]y[i])  2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i; 1320 } 1321 if (*chisq <= 0.0) *chisq = 1.0e308; // don't let chisq=0 through yfit being all ve 1322 break; 1323 1324 default: 1325 return 3; 1326 /* break; */ // (unreachable code) 1327 } 1328 1329 /* Fill in the symmetric side */ 1330 for (j=1; j<mfit; j++) 1331 for (k=0; k<j; k++) 1332 alpha[k][j] = alpha[j][k]; 1333 1334 return 0; 1335 } 1336 1337 1338 /* And this is the variant which handles an instrument response. */ 793 1339 /* We assume that the function values are sensible. */ 794 1340 … … 855 1401 if (paramfree[j]) mfit++; 856 1402 1403 //TODO ARG first change: 1404 // /* Initialise (symmetric) alpha, beta */ 1405 // for (j=0; j<mfit; j++) { 1406 // for (k=0; k<=j; k++) 1407 // alpha[j][k] = 0.0; 1408 // beta[j] = 0.0; 1409 // } 1410 857 1411 /* Calculation of the fitting data will depend upon the type of 858 1412 noise and the type of instrument response */ … … 931 1485 /* OK, now we've got our (possibly convolved) data, we can do the 932 1486 rest almost exactly as above. */ 1487 //TODO ARG this new section of code is just temporary; o'wise have to define all these variables at the top of the function 933 1488 { 934 float alpha_weight[ 256]; //TODO establish maximum # bins and use elsewhere (#define)935 float beta_weight[ 256];1489 float alpha_weight[MAXBINS]; 1490 float beta_weight[MAXBINS]; 936 1491 int q; 937 1492 float weight; … … 942 1497 float beta_sum; 943 1498 float dy_dparam_k_i; 944 float *alpha_weight_ptr;945 float *beta_weight_ptr;946 1499 947 1500 *chisq = 0.0f; … … 950 1503 case NOISE_CONST: 951 1504 { 952 float *alpha_ptr = &alpha_weight[fit_start];953 float *beta_ptr = &beta_weight[fit_start];954 1505 for (q = fit_start; q < fit_end; ++q) { 955 1506 (*pdy_dparam_conv)[q][0] = 1.0f; 956 1507 yfit[q] += param[0]; 957 1508 dy[q] = y[q]  yfit[q]; 958 weight = 1.0f; 959 *alpha_ptr++ = weight; // 960 //alpha_weight[q] = weight; // 1 1509 weight = 1.0f; //TODO ARG this should be 1.0f / sig[0] ! 1510 alpha_weight[q] = weight; // 1 961 1511 weight *= dy[q]; 962 //printf("beta weight %d %f is y %f  yfit %f\n", q, weight, y[q], yfit[q]); 963 *beta_ptr++ = weight; // 964 //beta_weight[q] = weight; // dy[q] 1512 beta_weight[q] = weight; // dy[q] 965 1513 weight *= dy[q]; 966 1514 *chisq += weight; // dy[q] * dy[q] … … 1068 1616 if (paramfree[j]) { 1069 1617 dot_product = 0.0f; 1070 alpha_weight_ptr = &alpha_weight[fit_start];1071 1618 if (0 == j_free) { 1072 beta_weight_ptr = &beta_weight[fit_start];1073 1619 // for all data 1074 1620 for (k = fit_start; k < fit_end; ++k) { 1075 1621 dy_dparam_k_i = (*pdy_dparam_conv)[k][i]; 1076 dot_product += dy_dparam_k_i * (*pdy_dparam_conv)[k][j] * (*alpha_weight_ptr++); //alpha_weight[k]; //TODO make it [i][k] and just *dy_dparam++ it.1077 beta_sum += dy_dparam_k_i * (*beta_weight_ptr++); ///beta_weight[k];1622 dot_product += dy_dparam_k_i * (*pdy_dparam_conv)[k][j] * alpha_weight[k]; //TODO make it [i][k] and just *dy_dparam++ it. 1623 beta_sum += dy_dparam_k_i * beta_weight[k]; 1078 1624 } 1079 1625 } … … 1081 1627 // for all data 1082 1628 for (k = fit_start; k < fit_end; ++k) { 1083 dot_product += (*pdy_dparam_conv)[k][i] * (*pdy_dparam_conv)[k][j] * (*alpha_weight_ptr++); //alpha_weight[k];1629 dot_product += (*pdy_dparam_conv)[k][i] * (*pdy_dparam_conv)[k][j] * alpha_weight[k]; 1084 1630 } 1085 1631 } // k loop … … 1098 1644 //else printf("param not free %d\n", i); 1099 1645 } // i loop 1100 1101 /* printf("i_free is %d j_free %d\n", i_free, j_free); 1102 //TODO try padding with zeroes; leftovers would affect AX=B solution 1103 for (i = i_free; i < nparam; ++i) { 1104 for (j = 0; j < nparam; ++j) { 1105 alpha[i][j] = alpha[j][i] = 0.0f; 1106 } 1107 beta[i] = 0.0f; 1108 }*/ 1109 1110 } 1646 } 1647 1111 1648 return 0; 1112 1649 } … … 1120 1657 noise models. They do not calculate alpha or beta. */ 1121 1658 1122 /* This is the variant which handles an instrument response. */ 1659 int GCI_marquardt_compute_fn_final(float x[], float y[], int ndata, 1660 noise_type noise, float sig[], 1661 float param[], int paramfree[], int nparam, 1662 void (*fitfunc)(float, float [], float *, float [], int), 1663 float yfit[], float dy[], float *chisq) 1664 { 1665 int i, j, mfit; 1666 float sig2i, dy_dparam[MAXFIT]; /* dy_dparam needed for fitfunc */ 1667 1668 for (j=0, mfit=0; j<nparam; j++) 1669 if (paramfree[j]) 1670 mfit++; 1671 1672 /* Calculation of the fitting data will depend upon the type of 1673 noise and the type of instrument response */ 1674 1675 /* There's no convolution involved in this function. This is then 1676 fairly straightforward, depending only upon the type of noise 1677 present. Since we calculate the standard deviation at every 1678 point in a different way depending upon the type of noise, we 1679 will place our switch(noise) around the entire loop. */ 1680 1681 switch (noise) { 1682 case NOISE_CONST: 1683 *chisq = 0.0; 1684 /* Summation loop over all data */ 1685 for (i=0; i<ndata; i++) { 1686 (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 1687 yfit[i] += param[0]; 1688 dy[i] = y[i]  yfit[i]; 1689 /* And find chi^2 */ 1690 *chisq += dy[i] * dy[i]; 1691 } 1692 1693 /* Now divide everything by sigma^2 */ 1694 sig2i = 1.0 / (sig[0] * sig[0]); 1695 *chisq *= sig2i; 1696 break; 1697 1698 case NOISE_GIVEN: /* This is essentially the NR version */ 1699 *chisq = 0.0; 1700 /* Summation loop over all data */ 1701 for (i=0; i<ndata; i++) { 1702 (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 1703 yfit[i] += param[0]; 1704 sig2i = 1.0 / (sig[i] * sig[i]); 1705 dy[i] = y[i]  yfit[i]; 1706 /* And find chi^2 */ 1707 *chisq += dy[i] * dy[i] * sig2i; 1708 } 1709 break; 1710 1711 case NOISE_POISSON_DATA: 1712 *chisq = 0.0; 1713 /* Summation loop over all data */ 1714 for (i=0; i<ndata; i++) { 1715 (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 1716 yfit[i] += param[0]; 1717 /* we still don't let the variance drop below 1 */ 1718 sig2i = (y[i] > 1 ? 1.0/y[i] : 1.0); 1719 dy[i] = y[i]  yfit[i]; 1720 /* And find chi^2 */ 1721 *chisq += dy[i] * dy[i] * sig2i; 1722 } 1723 break; 1724 1725 case NOISE_POISSON_FIT: 1726 *chisq = 0.0; 1727 // Summation loop over all data 1728 for (i=0; i<ndata; i++) { 1729 (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 1730 yfit[i] += param[0]; 1731 // we still don't let the variance drop below 1 1732 sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); 1733 dy[i] = y[i]  yfit[i]; 1734 // And find chi^2 1735 *chisq += dy[i] * dy[i] * sig2i; 1736 } 1737 break; 1738 1739 case NOISE_MLE: 1740 *chisq = 0.0; 1741 /* Summation loop over all data */ 1742 for (i=0; i<ndata; i++) { 1743 (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 1744 yfit[i] += param[0]; 1745 // dy[i] = y[i]  yfit[i]; 1746 1747 /* And find chi^2 */ 1748 // sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); 1749 // *chisq += dy[i] * dy[i] * sig2i; 1750 if (yfit[i]<=0.0) 1751 ; // do nothing 1752 else if (y[i]==0.0) 1753 *chisq += 2.0*yfit[i]; // to avoid NaN from log 1754 else 1755 *chisq += 2.0*(yfit[i]y[i])  2.0*y[i]*log(yfit[i]/y[i]); // was dy[i] * dy[i] * sig2i; 1756 } 1757 if (*chisq <= 0.0) *chisq = 1.0e308; // don't let chisq=0 through yfit being all ve 1758 break; 1759 1760 case NOISE_GAUSSIAN_FIT: 1761 *chisq = 0.0; 1762 // Summation loop over all data 1763 for (i=0; i<ndata; i++) { 1764 (*fitfunc)(x[i], param, &yfit[i], dy_dparam, nparam); 1765 yfit[i] += param[0]; 1766 sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); 1767 dy[i] = y[i]  yfit[i]; 1768 // And find chi^2 1769 *chisq += dy[i] * dy[i] * sig2i; 1770 } 1771 break; 1772 1773 default: 1774 return 3; 1775 /* break; */ // (unreachable code) 1776 } 1777 1778 return 0; 1779 } 1780 1781 1782 /* And this is the variant which handles an instrument response. */ 1123 1783 /* We assume that the function values are sensible. Note also that we 1124 1784 have to be careful about which points which include in the … … 1132 1792 float param[], int paramfree[], int nparam, 1133 1793 void (*fitfunc)(float, float [], float *, float [], int), 1134 float yfit[], float dy[], float *chisq, 1794 float yfit[], float dy[], float *chisq, 1135 1795 float **pfnvals, float ***pdy_dparam_pure, float ***pdy_dparam_conv, 1136 1796 int *pfnvals_len, int *pdy_dparam_nparam_size) … … 1187 1847 for (i=0; i<ndata; i++) { 1188 1848 int convpts; 1189 1849 1190 1850 /* We wish to find yfit = fnvals * instr, so explicitly: 1191 1851 yfit[i] = sum_{j=0}^i fnvals[ij].instr[j] … … 1222 1882 dy_dparam_conv[i], nparam); 1223 1883 } 1224 1884 1225 1885 /* OK, now we've got our (possibly convolved) data, we can do the 1226 1886 rest almost exactly as above. */ … … 1293 1953 case NOISE_POISSON_FIT: 1294 1954 *chisq = 0.0; 1295 // Summation loop over all data 1955 // Summation loop over all data 1296 1956 for (i=0; i<fit_start; i++) { 1297 1957 yfit[i] += param[0]; … … 1301 1961 yfit[i] += param[0]; 1302 1962 dy[i] = y[i]  yfit[i]; 1303 // And find chi^2 1304 // we still don't let the variance drop below 1 1963 // And find chi^2 1964 // we still don't let the variance drop below 1 1305 1965 sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); 1306 1966 *chisq += dy[i] * dy[i] * sig2i; … … 1311 1971 } 1312 1972 break; 1313 1973 1314 1974 case NOISE_MLE: // for the final chisq report a normal chisq measure for MLE 1315 1975 *chisq = 0.0; 1316 // Summation loop over all data 1976 // Summation loop over all data 1317 1977 for (i=0; i<fit_start; i++) { 1318 1978 yfit[i] += param[0]; … … 1322 1982 yfit[i] += param[0]; 1323 1983 dy[i] = y[i]  yfit[i]; 1324 // And find chi^2 1984 // And find chi^2 1325 1985 // sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); 1326 1986 if (yfit[i]<=0.0) … … 1340 2000 case NOISE_GAUSSIAN_FIT: 1341 2001 *chisq = 0.0; 1342 // Summation loop over all data 2002 // Summation loop over all data 1343 2003 for (i=0; i<fit_start; i++) { 1344 2004 yfit[i] += param[0]; … … 1348 2008 yfit[i] += param[0]; 1349 2009 dy[i] = y[i]  yfit[i]; 1350 // And find chi^2 2010 // And find chi^2 1351 2011 sig2i = (yfit[i] > 1 ? 1.0/yfit[i] : 1.0); 1352 2012 *chisq += dy[i] * dy[i] * sig2i; … … 1357 2017 } 1358 2018 break; 1359 2019 1360 2020 default: 1361 2021 return 3; … … 1373 2033 // was DoEcfFittingEngine() included in EcfSingle.c by PRB, 3.9.03 1374 2034 1375 int GCI_marquardt_fitting_engine(float xincr, float *trans, int ndata, int fit_start, int fit_end, 2035 int GCI_marquardt_fitting_engine(float xincr, float *trans, int ndata, int fit_start, int fit_end, 1376 2036 float prompt[], int nprompt, 1377 2037 noise_type noise, float sig[], … … 1388 2048 if (ecf_exportParams) ecf_ExportParams_OpenFile (); 1389 2049 1390 // All of the work is done by the ECF module 2050 // All of the work is done by the ECF module 1391 2051 ret = GCI_marquardt_instr(xincr, trans, ndata, fit_start, fit_end, 1392 2052 prompt, nprompt, noise, sig, … … 1394 2054 fitted, residuals, covar, alpha, &local_chisq, 1395 2055 chisq_delta, chisq_percent, erraxes); 1396 2056 1397 2057 // changed this for version 2, did a quick test with 2150ps_200ps_50cts_450cts.ics to see that the results are the same 1398 2058 // NB this is also in GCI_SPA_1D_marquardt_instr() and GCI_SPA_2D_marquardt_instr() 1399 2059 oldChisq = 3.0e38; 1400 while (local_chisq>chisq_target && (local_chisq<oldChisq) && tries< (MAXREFITS*3))2060 while (local_chisq>chisq_target && (local_chisq<oldChisq) && tries<MAXREFITS) 1401 2061 { 1402 oldChisq = local_chisq; 2062 oldChisq = local_chisq; 1403 2063 tries++; 1404 2064 ret += GCI_marquardt_instr(xincr, trans, ndata, fit_start, fit_end, … … 1407 2067 fitted, residuals, covar, alpha, &local_chisq, 1408 2068 chisq_delta, chisq_percent, erraxes); 1409 } 1410 /* if (local_chisq<=chisq_target) { 1411 printf("local_chisq %f target %f\n", local_chisq, chisq_target); 1412 } 1413 if (local_chisq>=oldChisq) { 1414 printf("local_chisq %f oldChisq %f\n", local_chisq, oldChisq); 1415 } 1416 if (tries>=(MAXREFITS*3)) { 1417 printf("tries is %d\n", tries); 1418 printf("local_chisq %f oldChisq %f\n", local_chisq, oldChisq); 1419 } 1420 1421 show_timing();*/ 2069 } 1422 2070 1423 2071 if (chisq!=NULL) *chisq = local_chisq; … … 1441 2089 } 1442 2090 2091 1443 2092 // Emacs settings: 1444 2093 // Local variables:
Note: See TracChangeset
for help on using the changeset viewer.