Forked from
mlz / Frida
2163 commits behind the upstream repository.
-
Wuttke, Joachim authoredWuttke, Joachim authored
manip.cpp 31.00 KiB
//**************************************************************************//
//* FRIDA: fast reliable interactive data analysis *//
//* manip.cpp: manipulate data (select, average, rearrange, ...) *//
//* (C) Joachim Wuttke 1990-, v2(C++) 2001- *//
//* http://frida.sourceforge.net *//
//**************************************************************************//
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <iostream>
#include <algorithm>
#include <boost/format.hpp>
#include <gsl/gsl_sort.h>
#include "mystd.h"
#include "olm.h"
#include "expr.h"
#include <ask_simple_value.h>
#include "manip.h"
#include "edif.h" // wg. ListZ()
#include "xax_lex.h"
#include "index.h"
using namespace std;
using boost::format;
// local routines:
namespace NManip {
void ScaRemoveConstantZ( POlo f );
int PtsExeRetEqui( const vector<double> x, uint nblock, double tol,
IndexSet *Sel);
}
//***************************************************************************//
//* restrict forthcoming operations *//
//***************************************************************************//
//! Restrict forthcoming operations
void NManip::ProtectSpecs( bool on )
{
NOlm::SelAssert();
string sel;
if( on )
sel = wask( "Protect which spectra from forthcoming operations" );
else
sel = "-";
NOlm::IterateO fiter;
POlo fio;
while( fio = fiter() ) {
fio->ProtectedSpecs = CList( sel, 0, fio->nJ()-1 );
fio->lDoc.push_back( "m/ " + sel );
}
}
//***************************************************************************//
//* manipulations on data points *//
//***************************************************************************//
//! Delete (or retain) points according to list of indices.
void NManip::PtsSelect( string del_or_ret )
{
NOlm::SelAssert();
string action;
bool sel_del;
if ( del_or_ret=="d" ){
action = "Delete";
sel_del = true;
} else if ( del_or_ret=="r" ){
action = "Retain";
sel_del = false;
} else
throw string( "FATAL: invalid del_or_ret" );
static string pSel = "";
pSel = wask( action + " which points", pSel );
CList pLis( pSel );
NOlm::IterateD fiter;
POld fin;
while( (fin = fiter()) ) {
POld fout( new COld( *fin ) );
fout->V.clear();
fout->lDoc.push_back( "mp" + del_or_ret + " " + pSel );
for ( uint j=0; j<fin->nJ(); j++ ) {
PSpec sin = fin->VS(j);
PSpec sout( new CSpec );
sout->copy_zentry( sin );
pLis.evaluate( 0, sin->size()-1 );
for ( uint i=0; i<sin->size(); i++ ) {
if ( sel_del ^ pLis.contains(i) ) {
if( sin->dy.size() )
sout->push_xyd( sin->x[i], sin->y[i], sin->dy[i] );
else
sout->push_xy( sin->x[i], sin->y[i] );
}
}
if ( sout->size() )
fout->V.push_back( sout );
}
if ( !fout->V.size() )
throw string( "no data point left" );
NOlm::OloAdd( fout, fiter.k() );
}
}
//! Bin points.
void NManip::PtsAvge()
{
NOlm::SelAssert();
static string groups;
groups = index_ask("Start bins at points", groups);
if (groups=="") return;
NOlm::IterateD fiter;
POld fin;
while( fin = fiter() ) {
POld fout( new COld( *fin ) );
fout->V.clear();
fout->lDoc.push_back("mpa " + groups);
for( uint j=0; j<fin->nJ(); j++ ) {
PSpec sin = fin->VS(j);
PSpec sout( new CSpec() );
sout->copy_zentry( sin );
IndexSet Groups;
if (Groups.parse(groups, 0, sin->size()))
throw "File " + strg(fiter.k()) + ", spectrum " + strg(j) + " invalid group";
if (!Groups.contains(0))
throw string( "Group must start at 0" );
Groups.insert(sin->size());
IndexSetIterator ISet;
ISet.reset(Groups);
uint i, igi, igf, ng;
while (ISet(&igi) && ISet(&igf)) {
// printf("DEBUG igi=%u igf=%u\n", igi, igf);
ISet.stepback();
double xg = 0, yg = 0, eg = 0;
for (i=igi; i<igf; ++i) {
xg += sin->x[i];
yg += sin->y[i];
if( sin->dy.size() )
eg += SQR( sin->dy[i] );
}
ng = igf - igi;
if( sin->dy.size() )
sout->push_xyd(xg/ng, yg/ng, sqrt(eg)/ng );
else
sout->push_xy(xg/ng, yg/ng );
}
fout->V.push_back(sout);
}
NOlm::OloAdd( fout, fiter.k() );
}
}
//! Sort spectrum according to given expression.
void NManip::PtsSort()
{
NOlm::SelAssert();
string expr = sask("Sort points according to");
if (expr=="") return;
PTree T;
user_xaxparse( expr.c_str(), &T );
NOlm::IterateD fiter;
POld fin;
while( fin = fiter() ) {
uint k, n, i;
k = fiter.k();
POld fout( new COld( *fin ) );
fout->V.clear();
fout->lDoc.push_back("mpo "+expr);
for ( uint j=0; j<fin->nJ(); j++ ) {
PSpec sin = fin->VS(j);
PSpec sout( new CSpec );
sout->copy_zentry( sin );
n = sin->size();
vector<double> v(n);
T->tree_vec_val( &v, k, j );
size_t P[n];
gsl_sort_index (P, &(v[0]), 1, n);
sout->resize( n, sin->dy.size() );
for (i=0; i<n; ++i) {
sout->x[i] = sin->x[P[i]];
sout->y[i] = sin->y[P[i]];
if( sin->dy.size() )
sout->dy[i] = sin->dy[P[i]];
}
fout->V.push_back( sout );
} NOlm::OloAdd( fout, fiter.k() );
}
}
//! Average points that have the same x.
void NManip::PtsAvgeEq(void)
{
NOlm::SelAssert();
COld fout;
NOlm::IterateD fiter;
POld fin;
while( fin = fiter() ) {
POld fout( new COld( *fin ) );
fout->V.clear();
fout->lDoc.push_back("mpq # average points when x equal");
for ( uint j=0; j<fin->nJ(); j++ ) {
const PSpec sin = fin->VS(j);
PSpec sout( new CSpec );
sout->copy_zentry( sin );
uint n = sin->size();
uint ii=0;
for ( uint i=0; i<n; i=ii ) {
double xg = sin->x[i];
double yg = sin->y[i];
double eg;
if( sin->dy.size() )
eg = SQR(sin->dy[i]);
uint ng = 1;
for ( ii=i+1; ii<n; ++ii ) {
if (sin->x[ii] != xg) break;
yg += sin->y[ii];
if( sin->dy.size() )
eg += SQR( sin->y[ii] );
++ng;
}
if( sin->dy.size() )
sout->push_xyd(xg, yg/ng, sqrt(eg)/ng );
else
sout->push_xy(xg, yg/ng);
}
fout->V.push_back( sout );
}
NOlm::OloAdd( fout );
}
}
//! Symmetrize spectrum.
void NManip::PtsSymmetrize()
{
NOlm::SelAssert();
static string expr = "0";
expr = sask( "Assume symmetry with mirror positioned at", expr );
if ( expr=="" ) return;
PTree T;
user_xaxparse( expr.c_str(), &T );
NOlm::IterateD fiter;
POld fin;
while( fin = fiter() ) {
uint k = fiter.k();
POld fout( new COld( *fin ) ); fout->V.clear();
fout->lDoc.push_back("mpsym "+expr);
uint n, i, il, ih, nl, nh;
for (uint j=0; j<fin->nJ(); j++) {
PSpec sin = fin->VS(j);
PSpec sout( new CSpec );
sout->copy_zentry( sin );
if( mystd::sorted( sin->x )!=1 )
throw string( "not sorted" );
double step;
if( !mystd::is_equidist( &step, sin->x ) )
throw string( "not equidistant" );
double centre;
T->tree_point_val( ¢re, k, j );
il = (int) ( (centre - sin->x[0])/step );
ih = il + 1;
n = sin->size();
nl = il + 1;
nh = n - ih;
if( nl<1 || nh<1 )
throw string( "nothing to symmetrize" );
for (i=0; i<(nl<nh?nl:nh); ++i) {
if( sin->dy.size() )
sout->push_xyd( ( -sin->x[il-i] + sin->x[ih+i] ) / 2,
( sin->y[il-i] + sin->y[ih+i] ) / 2,
sqrt( SQR(sin->dy[il-i]) +
SQR(sin->dy[ih+i]) ) / 2 );
else
sout->push_xy( ( -sin->x[il-i] + sin->x[ih+i] ) / 2,
( sin->y[il-i] + sin->y[ih+i] ) / 2 );
}
if( nl<nh ){
for (i=nl; i<nh; ++i)
if( sin->dy.size() )
sout->push_xyd( sin->x[ih+i], sin->y[ih+i],
sin->dy[ih+i] );
else
sout->push_xy( sin->x[ih+i], sin->y[ih+i] );
} else {
for (i=nh; i<nl; ++i)
if( sin->dy.size() )
sout->push_xyd( -sin->x[il-i], sin->y[il-i],
sin->dy[il-i] );
else
sout->push_xy( -sin->x[il-i], sin->y[il-i] );
}
sout->z = sin->z;
fout->V.push_back(sout);
}
NOlm::OloAdd( fout, fiter.k() );
}
}
//! Delete error bars.
void NManip::PtsErrorDelete()
{
NOlm::SelAssert();
NOlm::IterateD fiter;
POld fin;
while( fin = fiter() ) {
POld fout( new COld( *fin ) );
fout->V.clear();
fout->lDoc.push_back("mdy-");
for ( uint j=0; j<fin->nJ(); j++ ) { PSpec sout( new CSpec( *(fin->VS(j)) ) );
sout->dy.clear();
fout->V.push_back(sout);
}
NOlm::OloAdd( fout );
}
}
//***************************************************************************//
//* manipulations on spectra *//
//***************************************************************************//
//! Remove z variable that does not vary with j.
void NManip::ScaRemoveConstantZ( POlo fio )
{
NOlm::SelAssert();
for( int iz=fio->nZ()-1; iz>=0; --iz ){
uint nj = fio->nJ();
if( nj<1 )
throw string( "file contains no spectrum" );
double z = fio->z(0,iz);
bool z_is_const = true;
for( uint j=1; j<nj; ++j ){
if( fio->z(j,iz) != z ){
z_is_const = false;
}
}
if( z_is_const ){
printf( " constant z%i becomes r%zu\n", iz, fio->RPar.size() );
fio->RPar.push_back( CParam( fio->ZCo[iz], z ) );
fio->remove_z( iz );
}
}
}
//! Delete (or retain) spectra according to list of indices.
void NManip::ScaSelect( string del_or_ret )
{
NOlm::SelAssert();
string action;
bool sel_ret;
if ( del_or_ret=="d" ){
action = "Delete";
sel_ret = false;
} else if ( del_or_ret=="r" ){
action = "Retain";
sel_ret = true;
} else
throw string( "BUG: invalid del_or_ret" );
CList JSel, JSelSorted;
static string jSel = "";
NOlm::JSelAsk( action + " which spectra", &jSel, &JSel );
NOlm::IterateO fiter;
POlo fin;
while( fin = fiter() ) {
POlo fout( fin->new_olo() );
JSel.evaluate( 0, fin->nJ()-1 );
fout->lDoc.push_back( "ms" + del_or_ret + " " + JSel.str());
for( int j=0; j<fin->nJ(); ++j ) {
if( !sel_ret ^ JSel.contains( j ) ){
fout->V.push_back( fin->new_zentry( j ) );
} }
ScaRemoveConstantZ( fout );
NOlm::OloAdd( fout, fiter.k() );
}
}
//! Bin spectra.
void NManip::ScaAvge()
{
NOlm::SelAssert();
CList JSel;
static string jSel = "";
NOlm::JSelAsk( "Start groups at spectra", &jSel, &JSel );
NOlm::IterateO fiter;
POlo fin;
while( fin = fiter() ) {
POld fd = P2D( fin );
POlc fc = P2C( fin );
POlo fout( fin->new_olo() );
JSel.evaluate( 0, fin->nJ()-1 );
if (JSel.V.size()<1 || JSel.V[0]!=0)
throw string( "spectrum selection must contain 0" );
fout->lDoc.push_back("msa " + JSel.str());
for(uint iv=0; iv<JSel.size(); iv++) {
uint ji = JSel.V[iv];
uint jf = ( iv<JSel.size()-1 ? JSel.V[iv+1] : fin->nJ() );
uint nz = fin->nZ();
vector<double> zout( nz );
for( uint iz=0; iz<nz; ++iz ){
uint mz = 0;
double z = 0;
for ( uint jj=ji; jj<jf; jj++ ) {
z += fin->z(jj,iz);
++mz;
}
zout[iz] = z/mz;
}
if( fd ) {
PSpec sout( new CSpec() );
bool with_dy = fd->VS(0)->dy.size();
uint n = fd->nPts(ji);
sout->resize( n, with_dy );
for ( uint jj=ji+1; jj<jf; jj++){
if ( fd->nPts(jj) != fd->nPts(ji) )
throw "spectrum " + strg(jj) +
" has other size than spectrum " + strg(ji);
for( uint i=0; i<fd->nPts(jj); ++i ){
if ( fd->VS(jj)->x[i] != fd->VS(ji)->x[i] )
throw "spectrum " + strg(jj) + " has other x[" +
strg(i) + "] than spectrum 0";
sout->x[i] = fd->VS(ji)->x[i];
}
}
for ( uint i=0; i<n; i++ ){
uint mj = 0;
double ym = 0;
double vm = 0;
for ( uint jj=ji; jj<jf; jj++) {
ym += fd->VS(jj)->y[i];
if( with_dy )
vm += SQR(fd->VS(jj)->dy[i] );
++mj;
}
sout->y[i] = mj==0 ? 0 : ym/mj;
if( with_dy )
sout->dy[i] = mj==0 ? 0 : sqrt(vm)/mj; }
sout->z = zout;
fout->V.push_back( sout );
} else if ( fc ) {
PCurve cout( new CCurve );
uint np = fc->nPar();
vector<double> p(np);
for ( uint ip=0; ip<np; ip++ ){
uint mj = 0;
double pm = 0;
for ( uint jj=ji; jj<jf; jj++) {
pm += fc->VC(jj)->P[ip];
++mj;
}
p[ip] = mj==0 ? 0 : pm/mj;
}
cout->P = p;
cout->z = zout;
fout->V.push_back( cout );
}
}
ScaRemoveConstantZ( fout );
NOlm::OloAdd( fout, fiter.k() );
}
}
//! Merge spectra, concatenating points.
void NManip::ScaJoin()
{
NOlm::SelAssert();
CList JSel;
static string jSel = "";
NOlm::JSelAsk( "Start groups at spectra", &jSel, &JSel );
NOlm::IterateD fiter;
POld fin;
while( fin = fiter() ) {
POld fout( new COld( *fin ) );
fout->V.clear();
JSel.evaluate( 0, fin->nJ()-1 );
if (JSel.V.size()<1 || JSel.V[0]!=0)
throw string( "spectrum selection must contain 0" );
fout->lDoc.push_back("msj " + JSel.str());
cerr << "WARNING: z just taken from first spectrum of each group\n";
for ( uint iv=0; iv<JSel.size(); iv++ ) {
uint ji, jf;
ji = JSel.V[iv];
jf = ( iv<JSel.size()-1 ? JSel.V[iv+1] : fin->nJ() );
PSpec sout( new CSpec( *(fin->VS(ji)) ) );
bool with_dy = fin->VS(ji)->dy.size();
for ( uint jj=ji+1; jj<jf; jj++){
for( uint i=0; i<fin->VS(jj)->size(); ++i )
if ( with_dy )
sout->push_xyd( fin->VS(jj)->x[i], fin->VS(jj)->y[i],
fin->VS(jj)->dy[i]);
else
sout->push_xy( fin->VS(jj)->x[i], fin->VS(jj)->y[i] );
}
fout->V.push_back(sout);
}
ScaRemoveConstantZ( fout );
NOlm::OloAdd( fout, fiter.k() );
}
}
//! After each span, insert given number of copies.
void NManip::ScaSpawn()
{
NOlm::SelAssert();
static int njj = 2;
int njjIn = iask("output will have how many copies of input spectrum", njj);
if (njjIn<1) return;
njj = njjIn;
NOlm::IterateO fiter;
POlo fin;
while( fin = fiter() ) {
POlo fout( fin->new_olo() );
fout->lDoc.push_back( "ms* " + strg(njj) );
fout->ZCo.push_back(CCoord("no-in-spawn", ""));
for( uint jj=0; jj<njj; ++jj ){
for( uint j=0; j<fin->nJ(); j++ ){
PZentry eout = fin->new_zentry( j );
eout->z.push_back( jj );
fout->V.push_back( eout );
}
}
NOlm::OloAdd( fout, fiter.k() );
}
}
//! Exchange x and z.
void NManip::ScaExch()
{
NOlm::SelAssert();
NOlm::IterateD fiter;
POld fin;
while( fin = fiter() ) {
POld fout( new COld (*fin ) );
fout->V.clear();
uint izco = fin->ZCo.size() - 1;
if ( izco==(uint)-1 )
throw string( "no input z coordinate" );
fout->lDoc.push_back(
"msx # exchanging x=" + fin->xco.str() +
" with z" + strg(izco) + "=" + fin->ZCo[izco].str() );
fout->xco = fin->ZCo[izco];
fout->ZCo[izco] = fin->xco;
if ( fin->nJ()<=0 )
throw string( "no spectra in file" );
vector<double> zcommon, xcommon;
uint ji, jf;
for ( ji=0; ji<fin->nJ(); ) {
zcommon.clear();
for ( uint i=0; i<izco; ++i)
zcommon.push_back(fin->V[ji]->z[i]);
for ( jf=ji+1; jf<fin->nJ(); jf++ )
for ( uint i=0; i<izco; ++i)
if (fin->V[jf]->z[i] != zcommon[i])
goto end_group;
end_group:
// build common grid: xcommon.clear();
for ( uint j=ji; j<jf; ++j )
for ( uint i=0; i<fin->VS(j)->size(); ++i )
xcommon.push_back( fin->VS(j)->x[i] );
sort( xcommon.begin(), xcommon.end() );
mystd::unique(&xcommon, 1e-180, 1e-20);
printf("group of spectra %u .. %u has %zu different x\n",
ji, jf-1, xcommon.size());
uint ii = 0; // a guess
for ( uint i=0; i<xcommon.size(); ++i ) {
PSpec sout( new CSpec );
sout->z = zcommon;
sout->z.push_back(xcommon[i]);
for ( uint j=ji; j<jf; ++j ) {
PSpec sin = fin->VS( j );
bool with_dy = sin->dy.size();
if( sin->x[ii]!=xcommon[i] ){ // guess failed
for (ii=0; ii<sin->size(); ++ii)
if( sin->x[ii]==xcommon[i] )
break;
}
if( ii<sin->size() ){
if( with_dy )
sout->push_xyd( sin->z[izco], sin->y[ii],
sin->dy[ii] );
else
sout->push_xy( sin->z[izco], sin->y[ii] );
}
}
fout->V.push_back( sout );
++ii;
}
ji = jf;
}
NOlm::OloAdd( fout, fiter.k() );
}
}
//! Sort spectra according to expression.
void NManip::ScaSortByExpr()
{
NOlm::SelAssert();
string expr = sask("Sort spectra according to");
if (expr=="") return;
PTree T;
user_xaxparse( expr.c_str(), &T );
NOlm::IterateO fiter;
POlo fin;
while( fin = fiter() ) {
POlo fout( fin->new_olo() );
fout->lDoc.push_back("mso "+expr);
uint nj = fin->nJ();
vector<double> v(nj);
for (uint j=0; j<nj; j++)
T->tree_point_val( &(v[j]), fiter.k(), j );
size_t pos[nj];
gsl_sort_index (pos, &(v[0]), 1, nj);
for (uint j=0; j<nj; j++)
fout->V.push_back( fin->new_zentry( pos[j] ) );
NOlm::OloAdd( fout, fiter.k() );
}
}
//! Sort spectra according to full z vector.
void NManip::ScaSortByZ()
{
NOlm::SelAssert();
NOlm::IterateO fiter;
POlo fin;
while( fin = fiter() ) {
POlo fout( fin->new_olo() );
for (uint j=0; j<fin->nJ(); j++)
fout->V.push_back( fin->new_zentry( j ) );
sort( fout->V.begin(), fout->V.end(), CompareZ );
NOlm::OloAdd( fout, fiter.k() );
}
}
//! Change order of z coordinates.
void NManip::ZExchange()
{
NOlm::SelAssert();
int nzmin, nzmax = 0;
char mod;
string com;
int num1, num2;
NOlm::IterateO fiter;
POlo fin;
while( fin = fiter() ) {
nzmax = max( (int)nzmax, (int)fin->ZCo.size() );
}
if ( nzmax<=1 )
throw string( "nothing to exchange" );
else if( nzmax==2 ){
mod = 'r';
printf( "exchanging two z coordinates\n" );
com = "exchanged two z coordinates";
num1 = 1;
nzmin = 2;
} else {
string sel;
sel = wask( "rotate(r,rr,...), revert(m), exchange(<number>)" );
if ( sel[0]=='r' ){
mod = 'r';
num1 = 1;
while( sel.length() > num1 ){
if( sel[num1++]!='r' )
throw string( "invalid selection" );
}
com = "rotated z coordinates by " + strg(num1);
nzmin = 2;
} else if( sel=="m" ){
mod = 'm';
com = "reverted z coordinates";
nzmin = 2;
} else {
if( sscanf(sel.c_str(), "%d", &num1)!=1 )
throw string( "invalid selection" );
if( num1>=nzmax )
throw "no more than " + strg(nzmax) + " z coordinates";
num2 = iask( "exchange with coordinate" ); if( num2>=nzmax )
throw "no more than " + strg(nzmax) + " z coordinates";
else if( num1==num1 )
throw string( "nothing to exchange" );
mod = 'x';
com = "exchanged z coordinates "+strg(num1)+" and "+strg(num2);
nzmin = max( num1, num2 )+1;
}
}
fiter.reset();
while( fin = fiter() ) {
POlo fout( fin->new_olo() );
fout->lDoc.push_back( com );
for (uint j=0; j<fin->nJ(); j++)
fout->V.push_back( fin->new_zentry( j ) );
int nz = fout->ZCo.size();
if ( nz<nzmin ){
cerr << "WARNING nz<nzmin\n";
} else if( mod=='r' ){
for( int irot=0; irot<num1; ++irot ){
for( int iz=1; iz<nz; ++iz ){
fout->ZCo[iz] = fin->ZCo[iz-1];
}
fout->ZCo[0] = fin->ZCo[nz-1];
}
} else if( mod=='m' ){
for( int iz=0; iz<nz; ++iz ){
fout->ZCo[iz] = fin->ZCo[nz-1-iz];
}
} else { // mod=='x'
fout->ZCo[num1] = fin->ZCo[num2];
fout->ZCo[num2] = fin->ZCo[num1];
}
for ( uint j=0; j<fin->V.size(); j++ ) {
if ( nz<nzmin ){
// do nothing
} else if( mod=='r' ){
for( int irot=0; irot<num1; ++irot ){
for( int iz=1; iz<nz; ++iz )
fout->V[j]->z[iz] = fin->V[j]->z[iz-1];
fout->V[j]->z[0] = fin->V[j]->z[nz-1];
}
} else if( mod=='m' ){
for( int iz=0; iz<nz; ++iz ){
fout->V[j]->z[iz] = fin->V[j]->z[nz-1-iz];
}
} else { // mod=='x'
fout->V[j]->z[num1] = fin->V[j]->z[num2];
fout->V[j]->z[num2] = fin->V[j]->z[num1];
}
}
NOlm::OloAdd( fout, fiter.k() );
}
}
//! Delete z coordinate.
void NManip::ZDelete()
{
NOlm::SelAssert();
string sel;
NEdif::ListZ();
sel = wask("Delete z coordinates (-=quit)");
if (sel=="-" || sel=="q") return;
NOlm::IterateO fiter;
POlo fin;
while( fin = fiter() ) {
POlo fout( fin->new_olo() );
fout->lDoc.push_back("mz- "+sel);
for (uint j=0; j<fin->nJ(); j++)
fout->V.push_back( fin->new_zentry( j ) );
IndexSet ISel;
if ( ISel.parse(sel, 0, fin->ZCo.size()) )
throw ( "invalid selection for file " + strg( fiter.k() ) );
IndexSetIterator II;
II.setback(ISel);
uint iz;
while (II(&iz)) {
fout->ZCo.erase(fout->ZCo.begin()+iz);
for (uint j=0; j<fin->nJ(); j++)
fout->V[j]->z.erase(fout->V[j]->z.begin()+iz);
}
NOlm::OloAdd( fout, fiter.k() );
}
}
//! Break file into several files, according to z0.
void NManip::ScaBreak()
{
NOlm::SelAssert();
NOlm::IterateO fiter;
POlo fin;
while( fin = fiter() ) {
if ( !(fin->ZCo.size()) )
throw string( "no z coordinate" );
CCoord zco = fin->ZCo[0];
// intermediate file, on which output files will be based:
POlo ftmp( fin->new_olo() );
ftmp->V.clear();
ftmp->lDoc.push_back("spectra -> files, eliminating " + zco.str());
ftmp->ZCo.erase( ftmp->ZCo.begin() );
for ( uint j=0; j<fin->nJ(); ) {
POlo fout( ftmp->new_olo() );
double zval = ftmp->z(j,0);
fout->RPar.push_back( CParam( zco, zval ) );
do {
PZentry eout( fin->new_zentry(j) );
eout->z.erase( eout->z.begin() );
fout->V.push_back( eout );
} while ( ++j<fin->nJ() && fin->z(j,0)==zval );
NOlm::OloAdd( fout );
}
}
}
//***************************************************************************//
//* manipulations on files *//
//***************************************************************************//
//! Merge files, concatenating spectra.
void NManip::FilMerge( const string& opts )
{
NOlm::SelAssert();
/* TODO
COld fout;
CCoord zco;
vector<CParam> RPar;
NOlm::IterateD fiter;
POld fin;
// parse options:
bool add_zk = false;
if ( opts=="" )
;
else if ( opts=="+" )
add_zk = true;
// x,y coordinates:
CCoord co;
fiter.reset();
co = fiter()->xco;
while( (fin = fiter()) )
if ( fin->xco != co )
throw string( "different x coordinates" );
fout->xco = co;
fiter.reset();
co = fiter()->yco;
while ((fin = fiter()))
if (fin->yco != co)
throw string( "different y coordinates" );
fout->yco = co;
// z coordinates:
vector<CCoord> ZCo;
fiter.reset();
ZCo = fiter()->ZCo;
while ((fin = fiter())) {
if (fin->ZCo.size() != ZCo.size())
throw string( "different numer of z coordinates" );
for (uint i=0; i<ZCo.size(); ++i)
if (fin->ZCo[i] != ZCo[i])
throw "different z"+strg(i)+" coordinate";
}
fout->ZCo = ZCo;
if ( add_zk )
fout->ZCo.push_back(CCoord("number-of-input-file", ""));
// check RPar:
fiter.reset();
RPar = fiter()->RPar; // set RPar from first input file
vector<bool> RParDiff( RPar.size(), false );
while ((fin = fiter())) {
if (fin->RPar.size() != RPar.size())
throw string( "different numer of r parameters");
for (uint i=0; i<RPar.size(); ++i)
if (fin->RPar[i].Co != RPar[i].Co )
throw string( "different r parameter");
for (uint i=0; i<RPar.size(); ++i)
if (fin->RPar[i].val != RPar[i].val )
RParDiff[i] = true;
}
for (uint i=0; i<RPar.size(); ++i) {
if( RParDiff[i] )
fout->ZCo.push_back( RPar[i].Co );
else
fout->RPar.push_back( RPar[i] );
}
// file name:
string fnam;
fiter.reset();
fnam = fiter()->name; while ((fin = fiter())) {
uint n2=fin->name.size();
if (n2<fnam.size())
fnam.erase(fnam.begin()+n2, fnam.end());
for (uint i=0; i<fnam.size(); ++i)
if (fnam[i]!=(fin->name)[i])
fnam[i] = ' ';
}
for (string::iterator cp = fnam.begin(); cp != fnam.end(); ++cp)
if (*cp==' ') fnam.erase(cp--);
if (fnam == "") {
static uint callNo=0;
fnam = str( format( "merge%03u" ) % callNo++ );
}
fnam = wask("Joined file name", fnam);
if (fnam=="")
return;
fout->name = fnam;
// doc lines:
fout->lDoc.push_back(fnam + " is merger of:");
fiter.reset();
while ((fin = fiter())) {
fout->lDoc.push_back(" " + fin->name);
for (uint i=0; i<fin->lDoc.size(); ++i)
fout->lDoc.push_back(" " + fin->lDoc[i]);
}
// now merge the data:
fiter.reset();
while ((fin = fiter())) {
for (uint j=0; j<fin->nJ(); j++) {
PSpec S( new CSpec( *(fin->VS(j)) ) );
S->z.clear();
if ( add_zk )
S->z.push_back(fiter.index());
for (uint i=0; i<RPar.size(); ++i)
if( RParDiff[i] )
S->z.push_back( fin->RPar[i].val );
fout->V.push_back(S);
}
}
NOlm::OloAdd( fout );
*/
}
//! Merge files with same number of spectra, concatenating spectra with same j.
void NManip::FilMergePointwise()
{
NOlm::SelAssert();
POld fout; // TODO = *fin;
throw string("TODO FilMergePointwise probably broken after V rewrite");
NOlm::IterateD fiter;
POld fin;
while( fin = fiter() ) {
if ( fin->nJ() != fout->nJ() )
throw string( "different # spectra" );
if ( fin->xco != fout->xco || fin->yco != fout->yco )
throw string( "different x/y coordinates" );
for ( uint j=0; j<fout->nJ(); ++j ) {
const PSpec sin = fin->VS(j);
PSpec sout = fout->VS(j); if ( sin->z.size() != sout->z.size() )
throw string( "different # z entries" );
for ( uint i=0; i<sout->z.size(); ++i)
if ( sin->z[i] != sout->z[i] )
throw string( "different # z values" );
for ( uint i=0; i<sin->size(); ++i)
sout->push_xy(sin->x[i], sin->y[i]);
}
fout->lDoc.push_back("appended pointwise " + fin->name);
}
NOlm::OloAdd( fout );
}