further correction of msa

pub/CHANGELOG

 Release 2.3.3c:
 
 - Bug fix:
-  - msa: forbid breakpoints that exceed nJ
+  - msa: handle file-dependent nJ, handle mj
 
 Release 2.3.3b of 12jan16:
 

pub/lib/manip.cpp

@@ -476,31 +476,33 @@ void NManip::slices_rebin()
         // JSel.evaluate( 0, fin->nJ()-1 );
         if ( JSel.size()<1 || JSel[0]!=0 )
             throw S("spectrum selection must contain 0");
-        if ( JSel.back()>= fin->nJ() )
-            throw "Invalid breakpoint "+S(JSel.back())+": file "+S(fiter.k())+
-                " has only "+S(fin->nJ())+" slices";
+        vector<int> BreakPoints;
+        copy_if( JSel.begin(), JSel.end(), back_inserter(BreakPoints),
+                 [&](int j){ return j<fin->nJ(); } );
+        BreakPoints.push_back( fin->nJ() );
         fout->log_action("msa " + triv::indices_to_s(JSel));
-        for ( int iv=0; iv<JSel.size(); iv++ ) {
-            int ji = JSel[iv];
-            int jf = iv<JSel.size()-1 ? JSel[iv+1] : fin->nJ();
+        for ( int iv=0; iv<BreakPoints.size()-1; iv++ ) {
+            int ji = BreakPoints[iv];
+            int jf = BreakPoints[iv+1];
+            int mj = jf-ji;
+            if( mj<=0 )
+                throw S("Empty group");
             int nz = fin->nZ();
             vector<double> zout( nz );
             for ( int iz=0; iz<nz; ++iz ) {
-                int mz = 0;
                 double z = 0;
                 for ( int jj=ji; jj<jf; jj++ ) {
                     auto pz = PCAST<const CObjNum>(fin->z(jj,iz));
                     if( !pz )
                         throw "z"+S(iz)+" is not a number";
                     z += pz->to_r();
-                    ++mz;
                 }
-                zout[iz] = z/mz;
+                zout[iz] = z/mj;
             }
             if ( fd ) {
                 PSpec sout( new CSpec() );
                 bool in_with_dy = fd->VS(ji)->has_dy();
-                bool out_with_dy = in_with_dy || jf>ji;
+                bool out_with_dy = in_with_dy || mj>1;
                 int n = fd->nPts(ji);
                 sout->resize( n, out_with_dy );
                 // x grids must be equal:
@@ -519,12 +521,9 @@ void NManip::slices_rebin()
                 }
                 // bin y:
                 for ( int i=0; i<n; i++ ) {
-                    int mj = 0;
                     double ym = 0;
-                    for ( int jj=ji; jj<jf; jj++) {
+                    for ( int jj=ji; jj<jf; jj++)
                         ym += fd->VS(jj)->y[i];
-                        ++mj;
-                    }
                     if ( mj!=jf-ji || mj<1 )
                         throw S("BUG: inconsistent number of binned spectra");
                     ym /= mj;