diff --git a/pub/lib/fit.cpp b/pub/lib/fit.cpp
index 66ef2db241cb29e35432d02965304e03e44653fa..af966fad838e135d1b5eeac4f8b3a0b56bdc1186 100644
--- a/pub/lib/fit.cpp
+++ b/pub/lib/fit.cpp
@@ -179,18 +179,13 @@ void NCurveFit::fitOneSpec( POlc fc, ROld fd, int k, int j, const lm_control_str
// handle fixed parameters:
int np = fc->nP;
- int npfree = 0;
+ vector<double> Par;
for ( int ip=0; ip<np; ++ip )
if( fc->VC(j)->ParAttr[ip]!='x' )
- ++npfree;
+ Par.push_back( C->P[ip] );
+ int npfree = Par.size();
if ( npfree<=0 )
return;
- vector<double> Par(npfree);
- int ipf=0;
- for ( int ip=0; ip<np; ip++ ) {
- if ( fc->VC(j)->ParAttr[ip]!='x' )
- Par[ipf++] = C->P[ip];
- }
// Levenberg-Marquardt routine from library lmfit:
@@ -204,10 +199,9 @@ void NCurveFit::fitOneSpec( POlc fc, ROld fd, int k, int j, const lm_control_str
// Fit quality metrics Q:
- // Q: status code
C->fitOutcome = (double) status.outcome;
- // Q: coefficient of determination R^2 (for linear regression)
+ // coefficient of determination R^2 (for linear regression)
double ysum1 = 0, ysum2 = 0;
for( int i=0; i<nd; ++i ){
ysum1 += D->y[i];
@@ -216,7 +210,7 @@ void NCurveFit::fitOneSpec( POlc fc, ROld fd, int k, int j, const lm_control_str
double ydev2 = ysum2 - SQR( ysum1 )/nd; // = nd * <(y-<y>)^2>
C->fitR2 = ydev2!=0 ? 1 - SQR(status.fnorm) / ydev2 : -1;
- // Q: mean deviation (chi^2 if fvec is weighed with 1/dy)
+ // mean deviation (chi^2 if fvec is weighed with 1/dy)
int nfreedom = nd - npfree - 1;
if( nfreedom<1 ) nfreedom = 1;
C->fitChi2 = SQR(status.fnorm) / nfreedom;
@@ -228,8 +222,7 @@ void NCurveFit::fitOneSpec( POlc fc, ROld fd, int k, int j, const lm_control_str
printf( "%3d", j );
for ( int i=0; i<npfree; i++ )
printf(" %12g", Par[i]);
- printf( " (%3d->%9.7g)",
- status.nfev, C->fitChi2 );
+ printf( " (%3d->%9.7g)", status.nfev, C->fitChi2 );
printf( "%s\n", lm_shortmsg[ status.outcome ] );
}
}
@@ -244,7 +237,6 @@ typedef struct {
ROld fd;
POlc fc;
int k;
- int j;
double timeout;
} FitGloDatTyp;
@@ -338,38 +330,37 @@ void NCurveFit::fitGlobal( POlc fc, ROld fd, int k, const lm_control_struct& con
for( int j=0; j<fc->nJ(); ++j )
if ( fc->V[j]->frozen )
throw "Global fits currently incompatible with frozen spectra";
- for( int j=0; j<fc->nJ(); ++j ) {
- FitDatTyp data;
+ FitGloDatTyp data;
data.fc = fc;
data.fd = fd;
data.k = k;
- data.j = j;
- PSpec D = fd->VS(j);
- PCurve C = fc->VC(j);
- int nd = D->size();
- if ( !nd ){
- printf( "Empty data set in spectrum %i\n", j );
- return;
- }
// handle fixed parameters:
int np = fc->nP;
- int npfree = 0;
- for ( int ip=0; ip<np; ++ip )
- if( fc->VC(j)->ParAttr[ip]!='x' )
- ++npfree;
+ vector<double> Par;
+ for ( int ip=0; ip<np; ++ip ) {
+ for( int j=1; j<fc->nJ(); ++j )
+ if( (fc->VC(0)->ParAttr[ip]=='g') ^ (fc->VC(j)->ParAttr[ip]=='g') )
+ throw "Global attribute must currently be uniform for all spectra";
+ if ( fc->VC(0)->ParAttr[ip]=='g' ) {
+ Par.push_back( fc->VC(0)->P[ip] );
+ } else {
+ for( int j=0; j<fc->nJ(); ++j )
+ if ( fc->VC(0)->ParAttr[ip]=='u' )
+ Par.push_back( fc->VC(j)->P[ip] );
+ }
+ }
+ int npfree = Par.size();
if ( npfree<=0 )
return;
- vector<double> Par(npfree);
- int ipf=0;
- for ( int ip=0; ip<np; ip++ ) {
- if ( fc->VC(j)->ParAttr[ip]!='x' )
- Par[ipf++] = C->P[ip];
- }
-
+
+ int nd = 0;
+ for( int j=0; j<fc->nJ(); ++j )
+ nd += fd->VS(j)->size();
+
// Levenberg-Marquardt routine from library lmfit:
- data.timeout = time(NULL) + omp_get_num_threads()*maxtime;
+ data.timeout = time(NULL) + maxtime;
lm_status_struct status;
lmmin( npfree, &(Par[0]), nd, &data, fitGloEvaluate, &control, &status );
@@ -377,34 +368,34 @@ void NCurveFit::fitGlobal( POlc fc, ROld fd, int k, const lm_control_struct& con
if( status.outcome==11 )
throw "exception in fit function";
- // Fit quality metrics Q:
-
- // Q: status code
- C->fitOutcome = (double) status.outcome;
-
- // Q: coefficient of determination R^2 (for linear regression)
- double ysum1 = 0, ysum2 = 0;
- for( int i=0; i<nd; ++i ){
- ysum1 += D->y[i];
- ysum2 += SQR( D->y[i] );
- }
- double ydev2 = ysum2 - SQR( ysum1 )/nd; // = nd * <(y-<y>)^2>
- C->fitR2 = ydev2!=0 ? 1 - SQR(status.fnorm) / ydev2 : -1;
+ for( int j=0; j<fc->nJ(); ++j ) {
+ fc->VC(j)->fitOutcome = (double) status.outcome;
+
+ // coefficient of determination R^2 (for linear regression)
+ double ysum1 = 0, ysum2 = 0;
+ for( int i=0; i<nd; ++i ){
+ double y = fd->VS(j)->y[i];
+ ysum1 += y;
+ ysum2 += SQR( y );
+ }
+ double ydev2 = ysum2 - SQR( ysum1 )/nd; // = nd * <(y-<y>)^2>
+ fc->VC(j)->fitR2 = ydev2!=0 ? 1 - SQR(status.fnorm) / ydev2 : -1;
- // Q: mean deviation (chi^2 if fvec is weighed with 1/dy)
- int nfreedom = nd - npfree - 1;
- if( nfreedom<1 ) nfreedom = 1;
- C->fitChi2 = SQR(status.fnorm) / nfreedom;
-
+ // mean deviation (chi^2 if fvec is weighed with 1/dy)
+ int nfreedom = nd - npfree - 1;
+ if( nfreedom<1 ) nfreedom = 1;
+ fc->VC(j)->fitChi2 = SQR(status.fnorm) / nfreedom;
+ }
+
// print results:
- printf( "%3d", k );
- printf( "%3d", j );
- for ( int i=0; i<npfree; i++ )
- printf(" %12g", Par[i]);
- printf( " (%3d->%9.7g)",
- status.nfev, C->fitChi2 );
- printf( "%s\n", lm_shortmsg[ status.outcome ] );
- } //j
+ for( int j=0; j<fc->nJ(); ++j ) {
+ printf( "%3d", k );
+ printf( "%3d", j );
+ for ( int i=0; i<npfree; i++ )
+ printf(" %12g", Par[i]);
+ printf( " (%3d->%9.7g)", status.nfev, fc->VC(j)->fitChi2 );
+ printf( "%s\n", lm_shortmsg[ status.outcome ] );
+ }
}
//**************************************************************************************************