diff --git a/TODO b/TODO
index 3a1ade3e6d0dab03a72c0590b0c44b3e1071d308..7546f3c1c956b45374ef8f5615bd919504bfb84c 100644
--- a/TODO
+++ b/TODO
@@ -1,6 +1,8 @@
== Varia and Unsorted ==
-
+cp should show quality of last fit
+latest action on pX should be part of history
+resol with bg vs fit with bg+resol/conv
upgrade yaml-cpp from 0.3 to 0.5
== Building ==
diff --git a/pub/lib/fit.cpp b/pub/lib/fit.cpp
index d35eab483f7b4b3268756ff81af38f58df7cf074..da738de8798fbea66afd7ee1c180f9860d873b70 100644
--- a/pub/lib/fit.cpp
+++ b/pub/lib/fit.cpp
@@ -10,6 +10,7 @@
#include <iostream>
#include <string>
+#include <omp.h>
#include <lmmin.h>
#include <readplus/ask.hpp>
@@ -206,81 +207,84 @@ void NCurveFit::fit( bool _allow_slow_conv )
#pragma omp parallel for
for ( int j=0; j<fc->nJ(); j++ ) {
- if ( fc->V[j]->frozen )
- continue;
- FitDatTyp data;
- data.fc = fc;
- data.fd = fd;
- data.k = fiter.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 );
- continue;
- }
+ try {
+ if ( fc->V[j]->frozen )
+ continue;
+ FitDatTyp data;
+ data.fc = fc;
+ data.fd = fd;
+ data.k = fiter.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 );
+ continue;
+ }
- // handle fixed parameters:
- int npfree = 0;
- for ( int ip=0; ip<np; ++ip )
- if( !fc->VC(j)->fixed[ip] )
- ++npfree;
- if ( npfree<=0 )
- continue;
- vector<double> Par(npfree);
- int ipf=0;
- for ( int ip=0; ip<np; ip++ ) {
- if ( !fc->VC(j)->fixed[ip] )
- Par[ipf++] = C->P[ip];
- }
+ // handle fixed parameters:
+ int npfree = 0;
+ for ( int ip=0; ip<np; ++ip )
+ if( !fc->VC(j)->fixed[ip] )
+ ++npfree;
+ if ( npfree<=0 )
+ continue;
+ vector<double> Par(npfree);
+ int ipf=0;
+ for ( int ip=0; ip<np; ip++ ) {
+ if ( !fc->VC(j)->fixed[ip] )
+ Par[ipf++] = C->P[ip];
+ }
- if ( verbosity ) printf("\n");
-
- // Levenberg-Marquardt routine from library lmfit:
+ // Levenberg-Marquardt routine from library lmfit:
- data.timeout = time(NULL) + maxtime;
+ data.timeout = time(NULL) + omp_get_num_threads()*maxtime;
- lm_status_struct status;
- lmmin( npfree, &(Par[0]), nd, &data,
- fitEvaluate, &control, &status );
+ lm_status_struct status;
+ lmmin( npfree, &(Par[0]), nd, &data,
+ fitEvaluate, &control, &status );
- if( status.outcome==11 )
- throw "exception in fit function";
+ if( status.outcome==11 )
+ throw "exception in fit function";
- // Fit quality metrics Q:
+ // Fit quality metrics Q:
- // Q: status code
- C->Quality[0] = (double) status.outcome;
+ // Q: status code
+ C->Quality[0] = (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->Quality[1] = ydev2!=0 ? 1 - SQR(status.fnorm) / ydev2 : -1;
+ // 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->Quality[1] = 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->Quality[2] = SQR(status.fnorm) / nfreedom;
-
- #pragma omp critical
- {
- // print results:
- if( fiter.size()>1 )
- printf( "%3d", data.k );
- printf( "%3d", j );
- for ( int i=0; i<npfree; i++ )
- printf(" %12g", Par[i]);
- printf( " (%3d->%9.7g)",
- status.nfev, C->Quality[mFitMetric] );
- printf( "%s\n", lm_shortmsg[ status.outcome ] );
- }
- }
- }
- allow_slow_conv = true;
+ // Q: mean deviation (chi^2 if fvec is weighed with 1/dy)
+ int nfreedom = nd - npfree - 1;
+ if( nfreedom<1 ) nfreedom = 1;
+ C->Quality[2] = SQR(status.fnorm) / nfreedom;
+
+#pragma omp critical
+ {
+ // print results:
+ if( fiter.size()>1 )
+ printf( "%3d", data.k );
+ printf( "%3d", j );
+ for ( int i=0; i<npfree; i++ )
+ printf(" %12g", Par[i]);
+ printf( " (%3d->%9.7g)",
+ status.nfev, C->Quality[mFitMetric] );
+ printf( "%s\n", lm_shortmsg[ status.outcome ] );
+ }
+ } catch (...) {
+#pragma omp critical
+ printf( "Caught unexpected exception"
+ " while fitting spectrum %i\n", j );
+ } // try
+ } // j
+ } // k
}
diff --git a/pub/lib/func.cpp b/pub/lib/func.cpp
index e67598b5e045ebd5ac67f204efb574b2fef6d558..20ae940e129f30ceadc6c71a8a186aa25dd1cf62 100644
--- a/pub/lib/func.cpp
+++ b/pub/lib/func.cpp
@@ -297,6 +297,8 @@ double func_rotdiff( double w, double tau, double qb )
// Rotational diffusion spectrum
{
double ret = 0;
+ if( qb<=0 )
+ return 0;
for( int l=1; l<=12; ++l ){
ret += (2*l+1) * pow( gsl_sf_bessel_jl( l, qb ), 2 ) *
l*(l+1)*tau/PI/(SQR(w*tau)+SQR(l*(l+1)));