corrections; todo

TODO

 == 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 ==

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
 }
 

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)));