59 #include "ops_cmdline.h" 78 Interp(
char *ind,
char *weight );
79 int size() {
return ind.size(); }
94 FILE *iin = fopen( indf,
"r" );
95 while( fgets( buf, 8192, iin )!=NULL ) {
97 char *l1 = strdup(buf);
98 char *ptr = strtok( l1,
" " );
99 while((ptr=strtok( NULL,
" " ))) n++;
102 int *ip =
new int[n];
103 istringstream oss( buf );
104 for(
int i=0; i<n; i++ )
108 FILE *win = fopen( wf,
"r" );
109 for(
int i=0; i<ind.size(); i++ ) {
110 fgets( buf, 8192, win );
111 float *wp =
new float[n];
112 istringstream oss( buf );
113 for(
int j=0; j<n; j++ )
138 fprintf( stderr,
"Y must be an IGB file\n\n" );
142 gzFile in1 = (gzFile)(h1->
fileptr());
143 gzFile in2 = (gzFile)(h2->
fileptr());
144 FILE* out = (FILE *)(hout->
fileptr());
147 for(
int i=0; i<h1->
t(); i++ ) {
152 bool skip_second_torg =
overlap( h1, h2 );
153 for(
int i=0; i<h2->
t(); i++ ) {
156 skip_second_torg =
false;
166 fwrite( x,
sizeof(
data_t), nx, out );
167 fwrite( y,
sizeof(
data_t), ny, out );
176 static Parser hParser;
177 static bool init=
false;
180 hParser.DefineVar(
"X", x );
181 hParser.SetExpr(expr);
185 hParser.Eval(result, n);
187 catch (Parser::exception_type &e)
189 std::cerr << e.GetMsg() << std::endl;
215 hParser.DefineVar(
"X", &a );
216 hParser.DefineVar(
"Y", &b );
217 hParser.SetExpr(expr);
221 for(
int i=0; i<n; i++ ) {
224 b = x1[nc>0 ? i%nc : i];
225 y[i] = hParser.Eval();
227 catch (Parser::exception_type &e)
229 std::cerr << e.GetMsg() << std::endl;
241 int older= (t%2) ? 0 : slsz;
242 int old = older ? 0 : slsz;
245 for(
int i=0;i<slsz;i++){
246 data_t d = (x[i]-odat[old+i])/dt;
247 fwrite( &d,
sizeof(
data_t), 1, out );
250 for(
int i=0;i<slsz;i++){
251 data_t d = (x[i]-odat[older+i])/dt/2.;
252 fwrite( &d,
sizeof(
data_t), 1, out );
255 for(
int i=0;i<slsz;i++)
256 odat[older+i] = x[i];
263 for(
int i=0;i<slsz;i++){
265 data_t d = (x[i]-old[i])/dt;
266 fwrite( &d,
sizeof(
data_t), 1, out );
276 for(
int i=0;i<slsz;i++){
277 data_t datum = x[i]*a + b*y[i%nc];
278 fwrite( &datum,
sizeof(
data_t), 1, out );
286 for(
int i=0;i<slsz;i++){
287 data_t datum = x[i]*a + b*y[i];
288 fwrite( &datum,
sizeof(
data_t), 1, out );
296 for(
int c=0; c<nc; c++ ) {
298 #pragma omp parallel for reduction( +:s ) 299 for(
int i=c;i<slsz;i+=nc)
301 sumf[c] = s/double(slsz);
303 fwrite( &sumf,
sizeof(sumf[0]), nc, out );
310 #pragma omp parallel for 311 for(
int i=0; i<sum.size(); i++){
312 sum[i] += fac * x[i];
318 #pragma omp parallel for 319 for(
int i=0; i<sum.size(); i++){
320 sum[i] += fac * (x[i] - mean[i]) * (x[i] - mean[i]);
327 for(
int i=0;i<slsz;i++){
329 fwrite( &datum,
sizeof(
data_t), 1, out );
339 for(
int i=0;i<slsz;i++){
341 if (std::isinf(datum)) {
345 if (std::isnan(datum)){
349 fwrite( &datum,
sizeof(
data_t), 1, out );
351 if (INF_FLAG) std::cout <<
"Infinite number(s) detected (and set to zero)" << endl;
352 if (NAN_FLAG) std::cout <<
"NaN(s) detected (and set to zero)" << endl;
362 for(
int i=0;i<slsz;i++){
363 data_t datum = x[i]/y[i%nc];
364 if (std::isinf(datum)) {
368 if (std::isnan(datum)){
372 fwrite( &datum,
sizeof(datum), 1, out );
374 if (INF_FLAG) std::cout <<
"Infinite number(s) detected (and set to zero)" << endl;
375 if (NAN_FLAG) std::cout <<
"NaN(s) detected (and set to zero)" << endl;
381 for(
int i=0;i<slsz;i++){
383 fwrite( &datum,
sizeof(
data_t), 1, out );
390 for(
int i=0;i<slsz;i++){
391 data_t datum = x[i]*y[i%nc];
392 fwrite( &datum,
sizeof(datum), 1, out );
400 for(
int i=0; i<nc; i++ ) min[i]=x[i];
401 #pragma omp parallel for 402 for(
int i=nc;i<slsz*nc;i++){
406 fwrite( min,
sizeof(
data_t), nc, out );
413 #pragma omp parallel for 414 for(
int i=0; i<n; i++ ){
425 #pragma omp parallel for 426 for(
int i=0; i<n; i++ ){
436 #pragma omp parallel for 437 for(
int i=0; i<n; i++ ){
448 for(
int i=0; i<nc; i++ ) max[i]=x[i];
449 #pragma omp parallel for 450 for(
int i=nc;i<slsz*nc;i++){
454 fwrite( max,
sizeof(
data_t), nc, out );
460 for(
int i=0;i<slsz;i++){
462 fwrite( &datum,
sizeof(
data_t), 1, out );
469 void Xbox(
data_t *dat,
int x,
int y, FILE *out,
int d,
int nc )
476 for(
int j=0; j<y; j++ )
477 for(
int i=0; i<x; i++ ) {
479 if( j<fm || j>y-fp-1 || i<fm || i>x-fp-1 )
482 for(
int m=j-fm; m<=j+fp; m++ )
483 for(
int n=i-fm; n<=i+fp; n++ )
487 fwrite( &datum,
sizeof(
data_t), 1, out );
494 int tfinal, FILE *out )
496 unsigned int curr = (t%3)*sz;
497 unsigned int old = (curr+2*sz)%(3*sz);
498 unsigned int older = (curr+sz)%(3*sz);
500 memcpy( scr.data()+curr, dat,
sizeof(
data_t)*sz );
503 fwrite( dat,
sizeof(
data_t), sz, out );
506 for(
int i=0; i<sz; i++ )
507 if( scr[old+i]>upper || scr[old+i]<lower ) {
508 if( scr[curr+i]>upper || scr[curr+i]<lower )
509 scr[old+i] = scr[older+i];
511 scr[old+i] = (scr[older+i]+scr[curr+i])/2.;
513 fwrite( scr.data()+old,
sizeof(
data_t), sz, out );
517 fwrite( dat,
sizeof(
data_t), sz, out );
531 for(
int i=0; i<n; i+=nc){
533 for(
int j=0; j<nc; j++ )
534 mag += x[i+j]*x[i+j];
536 fwrite( &mag,
sizeof(mag), 1, out );
553 for(
int i=0; i<n; i+=nc){
555 for(
int j=0; j<nc; j++ )
556 dot += x[i+j]*y[k*i+j];
557 fwrite( &dot,
sizeof(dot), 1, out );
581 else if( s==1 || s==2 )
584 #pragma omp parallel for 585 for(
int i=0; i<n; i++ )
586 runInt[i] += x[i]*scale;
596 #pragma omp parallel for 597 for(
int i=0; i<n; i++ )
599 reduce[i] += dt*x[i]/2.;
601 reduce[i] += dt*x[i];
618 for(
int i=0; i<idx.size(); i++ ) {
620 for(
int j=0; j<n; j++ )
621 datum += dat[idx[i][j]]*w[i][j];
622 fwrite( &datum,
sizeof(
data_t), 1, out );
641 vector<data_t> slice( slice_mem );
643 for(
int t=0; t< h->
t(); t++ ) {
645 for(
int j=0; j<slice_mem; j++ ){
646 data_t datum = (slice[j]-min[j])/(max[j]-min[j]);
647 fwrite( &datum,
sizeof(
data_t), 1, out );
672 data_t *extrem = big+2*sz;
674 if( !t ) memcpy( old, dat,
sizeof(
data_t)*sz );
676 #pragma omp parallel for 677 for(
int i=0; i<sz; i++ ) {
678 data_t dV = dat[i] - old[i];
679 if( dV*change[i] < 0 ) {
680 data_t pp = fabs( extrem[i]-old[i] );
681 if( pp > maxpp[i] ) maxpp[i] = pp;
682 change[i] = copysign( 1., dV );
684 }
else if( change[i]==0. && dV != 0 ) {
685 change[i] = copysign( 1., dV );
688 memcpy( old, dat,
sizeof(
data_t)*sz );
696 for(
int i=0; i<n; i++ )
698 fwrite( s+i,
sizeof(
data_t), 1, out );
700 fprintf( out,
"%g\n", s[i] );
710 fwrite( &datum,
sizeof(
data_t), 1, out );
712 fprintf( out,
"%g\n", i );
718 void ravgX_ta(
data_t *dat1,
int sz,
int t,
int a, vector<data_t>& tmp, FILE *out )
720 vector <data_t>
sum( sz, 0 );
722 #pragma omp parallel for 723 for(
int i=0; i<sz; i++ ) {
725 tmp[i+(t%a)*sz] = dat1[i];
727 for(
int j=0; j<a; j++ ) {
728 sum[i] += tmp[i+j*sz];
747 while( gzgets( in, buf, 1024 ) != Z_NULL &&
748 sscanf( buf,
"%f %f %f %f", datum, datum+1, datum+2, datum+3 )==nc )
749 for(
int i=0; i<nc; i++ ) y.push_back( datum[i] );
762 for(
int i=0; i<n; i++ ) {
774 if( !getenv(
"OMP_NUM_THREADS") )
775 omp_set_num_threads(1);
778 gengetopt_args_info opts;
783 if (cmdline_parser (argc, argv, &opts) != 0)
787 gzFile in1 = (!opts.inputs_num || !strcmp(opts.inputs[0],
"-")) ?
788 gzdopen( dup(0),
"r" ):
789 gzopen( opts.inputs[0],
"r" );
791 cerr <<
"File not found: " << opts.inputs[0] << endl;
796 if( h1->
read(
true) & 0xff ) {
797 cerr <<
"File 1 is not a proper IGB file: "<< opts.inputs[0] << endl;
802 if( opts.interpolate_flag )
803 interp =
new Interp( opts.inputs[1], opts.inputs[2] );
808 switch( opts.op_arg ) {
809 case op_arg_aXPLUS_bY:
812 case op_arg_joinXY_t:
813 case op_arg_joinXY_i:
818 else if( opts.expression_given )
819 if( strchr( opts.expression_arg,
'Y' ) )
830 if( opts.inputs_num<2 ) {
831 cerr <<
"A 2nd file is needed : " << endl;
834 if( (in2=gzopen( opts.inputs[1],
"r" ))==NULL ) {
835 cerr <<
"file 2 not found : " << opts.inputs[1] << endl;
840 if( h2->
read(
true) & 0xff ) {
849 else if( opts.op_arg != op_arg_joinXY_i ) {
850 cerr <<
"incompatile file sizes! Y must have " << h1->
t()*h1->
num_components()
852 <<
" entries" << endl;
858 cerr <<
"Number of data components do not match" << endl;
865 if( h2->
t()!=1 && opts.op_arg != op_arg_joinXY_t ) {
866 cerr <<
"Second IGB file time does not match" << endl;
870 }
else if( h2->
slice_sz()==1 && h2->
t()==h1->
t() )
872 else if ( h2->
t()==h1->
t() )
874 else if ( h2->
slice_sz()==1 && h2->
t()==1 )
876 else if( opts.op_arg != op_arg_joinXY_i ) {
877 cerr <<
"incompatile IGB file sizes!" << endl;
886 string fname = opts.output_file_arg;
890 if( fname.length()<5 || fname.substr(fname.length()-4,4)!=
".igb" )
891 fname += opts.extension_arg;
892 out = fopen( fname.c_str(),
"w" );
897 switch( opts.op_arg ) {
898 case op_arg_ravgX_ta:
899 hout->t(h1->
t()-opts.a_arg+1);
901 case op_arg_joinXY_i:
909 case op_arg_joinXY_t:
910 hout->type(h1->
type());
940 if( opts.interpolate_flag ) {
941 hout->x( interp->
size() );
946 hout->fileptr( out );
949 if( opts.op_arg == op_arg_joinXY_t ) {
956 std::vector<std::vector<float>> dat3;
962 std::vector<data_t> t_reduce;
963 std::vector<data_t> big_scratch;
964 switch( opts.op_arg ) {
971 t_reduce.assign( h1_slice_comps, 0 );
974 t_reduce.assign( h1_slice_comps, 0 );
975 big_scratch.assign( h1_slice_comps, 0 );
977 case op_arg_ravgX_ta :
978 big_scratch.assign( h1_slice_comps*opts.a_arg, 0 );
980 case op_arg_sclip_ab :
981 big_scratch.assign( h1_slice_comps*3, 0 );
984 big_scratch.assign( h1_slice_comps, 0 );
988 big_scratch.assign( h1_slice_comps*3, 0 );
991 big_scratch.assign( h1_slice_comps*3, 0 );
992 t_reduce.assign( h1_slice_comps, 0 );
995 if( opts.expression_given )
996 big_scratch.assign( h1_slice_comps*opts.a_arg, 0 );
999 for(
int t=0; t<h1->
t(); t++ ) {
1005 if(opts.op_arg == op_arg_stdD_i)
1006 dat3.push_back(dat1);
1008 if( Yigb && (Ytype&
Tbit || !t) )
1011 float *Yptr = Y.data();
1017 if( opts.op_arg==op_arg_minX_i || opts.op_arg==op_arg_maxX_i ) {
1019 }
else if( opts.op_arg==op_arg_norm_i ) {
1020 copy( dat1.begin(), dat1.end(), big_scratch.begin() );
1021 copy( dat1.begin(), dat1.end(), big_scratch.begin()+h1_slice_comps );
1022 }
else if( opts.op_arg==op_arg_t_maxX) {
1027 if( opts.expression_given ){
1029 eval_expr( dat1.data(), h1_slice_comps, opts.expression_arg, out, big_scratch.data() );
1032 eval_expr( dat1.data(), Yptr, h1_slice_comps, opts.expression_arg,
1037 if( opts.interpolate_given )
1040 switch( opts.op_arg ) {
1041 case op_arg_aXPLUS_bY:
1043 aXPLUS_bY(dat1.data(), Yptr, opts.a_arg, opts.b_arg,
1044 h1_slice_comps,out);
1046 aXPLUS_bY(dat1.data(), Yptr, opts.a_arg, opts.b_arg,
1049 case op_arg_aXPLUS_b:
1050 aXPLUS_b(dat1.data(),opts.a_arg,opts.b_arg,h1_slice_comps,out);
1054 XdivY(dat1.data(), Yptr, h1_slice_comps, out);
1060 XmulY( dat1.data(), Yptr, h1_slice_comps, out);
1074 minXi( dat1.data(), t_reduce.data(), t_reduce.size() );
1077 maxXi( dat1.data(), t_reduce.data(), t_reduce.size() );
1080 t_maxX( dat1.data(), big_scratch.data(), t_reduce.data(), t_reduce.size(), t );
1084 avgXi( dat1.data(), t_reduce, 1.0f);
1086 case op_arg_ravgX_ta :
1087 ravgX_ta( dat1.data(), h1_slice_comps, t, opts.a_arg, big_scratch, out );
1092 case op_arg_sclip_ab :
1093 sclip_ab( dat1.data(), h1_slice_comps, opts.a_arg, opts.b_arg, big_scratch, t, h1->
t(), out );
1095 case op_arg_diff_f :
1098 case op_arg_diff_c :
1099 center_diff( dat1.data(), h1_slice_comps, big_scratch, t, h1->
inc_t(), out );
1101 case op_arg_joinXY_i:
1103 join_space( dat1.data(), h1_slice_comps, Yptr,
1106 join_space( dat1.data(), h1_slice_comps, Y.data(), Y.size(), out);
1115 case op_arg_intX_i :
1116 case op_arg_intX_t :
1118 int lastT = h1->
t()-1;
1119 int lastSimp = lastT - (lastT%3);
1120 if( t<=lastSimp && lastSimp!=0 )
1121 simp38( dat1.data(), h1_slice_comps, t==0||t==lastSimp?-1:t, t_reduce.data(),
1123 if( t>=lastSimp && lastSimp!=lastT )
1124 midpt( dat1.data(), h1_slice_comps, t==lastSimp||t==lastT, t_reduce.data(),
1127 if( opts.op_arg == op_arg_intX_t )
1128 fwrite( t_reduce.data(),
sizeof(t_reduce[0]), t_reduce.size(), out );
1131 minXi( dat1.data(), big_scratch.data(), h1_slice_comps );
1132 maxXi( dat1.data(), big_scratch.data()+h1_slice_comps, h1_slice_comps );
1135 maxP2P( dat1.data(), big_scratch.data(), t_reduce.data(), h1_slice_comps, t );
1143 switch( opts.op_arg ) {
1146 t_reduce[i] /= (
data_t)(h1->
t());
1151 t_reduce[i] /= (
data_t)(h1->
t());
1152 std::vector<data_t> mean;
1153 mean.assign(t_reduce.begin(), t_reduce.end());
1154 std::fill(t_reduce.begin(), t_reduce.end(), 0);
1155 for(
int t=0; t<h1->
t(); t++ ) {
1157 varXi( dat3[t].data(), mean.data(), t_reduce, 1.0f);
1160 t_reduce[i] = std::sqrt(t_reduce[i] / ((
data_t)(h1->
t() - 1)));
1164 normalize( h1, out, big_scratch.data(), big_scratch.data()+h1_slice_comps );
1167 for(
int i=0; i<t_reduce.size(); i++ )
1168 t_reduce[i] = h1->
org_t() + t_reduce[i]*h1->
inc_t();
1172 if( t_reduce.size() )
constexpr T min(T a, T b)
void read_vector(vector< data_t > &y, gzFile in, int nc)
void XdivY(data_t *x, data_t *y, int slsz, FILE *out)
void midpt(T *x, int n, bool end, data_t *reduce, float dt)
void normalize(IGBheader *h, FILE *out, T *min, T *max)
normalize pixels over all time
void aXPLUS_bY(data_t *x, data_t *y, data_t a, data_t b, int slsz, int nc, FILE *out)
void magnitude(data_t *x, int n, int nc, FILE *out)
determine magnitude at each point
void eval_expr(data_t *x, int n, char *expr, FILE *out, data_t *result)
constexpr T max(T a, T b)
void center_diff(data_t *x, int slsz, vector< data_t > &odat, int t, float dt, FILE *out)
void aXPLUS_b(data_t *x, data_t a, data_t b, int slsz, FILE *out)
void interpolate(data_t *dat, vector< int *> &idx, vector< float *> &w, int n, FILE *out)
interpolate within points
void forward_diff(data_t *x, int slsz, vector< data_t > &old, int t, float dt, FILE *out)
void avgXi(data_t *x, vector< data_t > &sum, const data_t fac)
void maxXi(data_t *x, data_t *mx, size_t n)
void output_scratch(data_t *s, FILE *out, int n, bool bin)
void XPLUS_a(data_t *x, data_t a, int slsz, FILE *out)
void joinXY(IGBheader *h1, IGBheader *h2, IGBheader *hout)
void XmulY(data_t *x, data_t *y, int slsz, FILE *out)
int main(int argc, char *argv[])
void ravgX_ta(data_t *dat1, int sz, int t, int a, vector< data_t > &tmp, FILE *out)
void maxXt(data_t *x, int slsz, int nc, float t, FILE *out)
std::string cmdline_str(int n, char *w[])
T sum(const vector< T > &vec)
Compute sum of a vector's entries.
void join_space(data_t *x, int nx, data_t *y, int ny, FILE *out)
Interp(char *ind, char *weight)
constructor
void varXi(data_t *x, data_t *mean, vector< data_t > &sum, const data_t fac)
bool overlap(IGBheader *first, IGBheader *second)
void maxP2P(data_t *dat, data_t *big, data_t *maxpp, int sz, int t)
find the maximum peak-to-peak value at each point
void simp38(T *x, int n, int t, data_t *runInt, float dt)
Simpson's 3/8 rule.
void t_maxX(data_t *x, data_t *mx, data_t *tm, size_t n, float t)
void dot(data_t *x, data_t *y, int n, int nc, FILE *out, bool pt)
perform dot product
void Xbox(data_t *dat, int x, int y, FILE *out, int d, int nc)
void avgXt(data_t *x, int slsz, int nc, float t, FILE *out)
void sclip_ab(data_t *dat, int sz, data_t upper, data_t lower, vector< data_t > &scr, int t, int tfinal, FILE *out)
void minXi(data_t *x, data_t *mn, size_t n)
double mag(const Point &vect)
vector magnitude
void minXt(data_t *x, int slsz, int nc, float t, FILE *out)