source: trunk/src/STIdxIter.cpp@ 2859

Last change on this file since 2859 was 2580, checked in by ShinnosukeKawakami, 12 years ago

hpc33 merged asap-trunk

File size: 12.8 KB
Line 
1#include <iostream>
2#include <casa/Utilities/GenSort.h>
3#include <casa/Arrays/Vector.h>
4#include <casa/Arrays/Matrix.h>
5#include <casa/Arrays/ArrayMath.h>
6#include <casa/Arrays/ArrayIO.h>
7#include <tables/Tables/ScalarColumn.h>
8#include "STIdxIter.h"
9
10namespace asap {
11
12IndexIterator::IndexIterator( IPosition &shape )
13 : niter_m( 0 )
14{
15 nfield_m = shape.nelements() ;
16 prod_m.resize( nfield_m ) ;
17 idx_m.resize( nfield_m ) ;
18 prod_m[nfield_m-1] = 1 ;
19 idx_m[nfield_m-1] = 0 ;
20 for ( Int i = nfield_m-2 ; i >= 0 ; i-- ) {
21 prod_m[i] = shape[i+1] * prod_m[i+1] ;
22 idx_m[i] = 0 ;
23 }
24 maxiter_m = prod_m[0] * shape[0] ;
25// cout << "maxiter_m=" << maxiter_m << endl ;
26}
27
28Bool IndexIterator::pastEnd()
29{
30 return niter_m >= maxiter_m ;
31}
32
33void IndexIterator::next()
34{
35 niter_m++ ;
36 uInt x = niter_m ;
37 for ( Int i = 0 ; i < nfield_m ; i++ ) {
38 idx_m[i] = x / prod_m[i] ;
39 //cout << "i=" << i << ": prod=" << prod_m[i]
40 // << ", idx=" << idx_m[i] << endl ;
41 x %= prod_m[i] ;
42 }
43}
44
45ArrayIndexIterator::ArrayIndexIterator( Matrix<uInt> &arr,
46 vector< vector<uInt> > idlist )
47 : arr_m( arr ),
48 pos_m( 1 )
49{
50 ncol_m = arr_m.ncolumn() ;
51 nrow_m = arr_m.nrow() ;
52 vector< vector<uInt> > l = idlist ;
53 if ( l.size() != ncol_m ) {
54 l.resize( ncol_m ) ;
55 for ( uInt i = 0 ; i < ncol_m ; i++ ) {
56 Vector<uInt> a( arr_m.column( i ).copy() ) ;
57 uInt n = genSort( a, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates ) ;
58 a.resize(n,True) ;
59 a.tovector( l[i] ) ;
60 }
61 }
62 idxlist_m = l ;
63 IPosition shape( ncol_m ) ;
64 for ( uInt i = 0 ; i < ncol_m ; i++ ) {
65 shape[i] = idxlist_m[i].size() ;
66 }
67// cout << "shape=" << shape << endl ;
68 iter_m = new IndexIterator( shape ) ;
69 current_m.resize( ncol_m ) ;
70}
71
72ArrayIndexIterator::~ArrayIndexIterator()
73{
74 delete iter_m ;
75}
76
77Vector<uInt> ArrayIndexIterator::current()
78{
79 Block<uInt> idx = iter_m->current() ;
80 for ( uInt i = 0 ; i < ncol_m ; i++ )
81 current_m[i] = idxlist_m[i][idx[i]] ;
82 return current_m ;
83}
84
85Bool ArrayIndexIterator::pastEnd()
86{
87 return iter_m->pastEnd() ;
88}
89
90ArrayIndexIteratorNormal::ArrayIndexIteratorNormal( Matrix<uInt> &arr,
91 vector< vector<uInt> > idlist )
92 : ArrayIndexIterator( arr, idlist )
93{
94 storage_m.resize( nrow_m ) ;
95}
96
97void ArrayIndexIteratorNormal::next()
98{
99 iter_m->next() ;
100}
101
102Vector<uInt> ArrayIndexIteratorNormal::getRows( StorageInitPolicy policy )
103{
104 Vector<uInt> v = current() ;
105 uInt len = 0 ;
106 uInt *p = storage_m.storage() ;
107 for ( uInt i = 0 ; i < nrow_m ; i++ ) {
108 if ( allEQ( v, arr_m.row( i ) ) ) {
109 *p = i ;
110 len++ ;
111 p++ ;
112 }
113 }
114 pos_m[0] = len ;
115 p = storage_m.storage() ;
116 return Vector<uInt>( pos_m, p, policy ) ;
117}
118
119ArrayIndexIteratorAcc::ArrayIndexIteratorAcc( Matrix<uInt> &arr,
120 vector< vector<uInt> > idlist )
121 : ArrayIndexIterator( arr, idlist )
122{
123 // storage_m layout
124 // length: ncol_m * (nrow_m + 1)
125 // 0~nrow_m-1: constant temporary storage that indicates whole rows
126 // nrow_m~2*nrow_m-1: temporary storage for column 0
127 // 2*nrow_m~3*nrow_m-1: temporary storage for column 1
128 // ...
129 storage_m.resize( arr_m.nelements()+nrow_m ) ;
130 // initialize constant temporary storage
131 uInt *p = storage_m.storage() ;
132 for ( uInt i = 0 ; i < nrow_m ; i++ ) {
133 *p = i ;
134 p++ ;
135 }
136 // len_m layout
137 // len[0]: length of temporary storage that incidates whole rows (constant)
138 // len[1]: number of matched row for column 0 selection
139 // len[2]: number of matched row for column 1 selection
140 // ...
141 len_m.resize( ncol_m+1 ) ;
142 p = len_m.storage() ;
143 for ( uInt i = 0 ; i < ncol_m+1 ; i++ ) {
144 *p = nrow_m ;
145 p++ ;
146// cout << "len[" << i << "]=" << len_m[i] << endl ;
147 }
148 // skip_m layout
149 // skip_m[0]: True if column 0 is filled by unique value
150 // skip_m[1]: True if column 1 is filled by unique value
151 // ...
152 skip_m.resize( ncol_m ) ;
153 for ( uInt i = 0 ; i < ncol_m ; i++ ) {
154 skip_m[i] = (Bool)(idxlist_m[i].size()==1) ;
155// cout << "skip_m[" << i << "]=" << skip_m[i] << endl ;
156 }
157 prev_m = iter_m->current() ;
158 p = prev_m.storage() ;
159 for ( uInt i = 0 ; i < ncol_m ; i++ ) {
160 *p = *p - 1 ;
161 p++ ;
162 }
163// cout << "prev_m=" << Vector<uInt>(IPosition(1,ncol_m),prev_m.storage()) << endl ;
164}
165
166void ArrayIndexIteratorAcc::next()
167{
168 prev_m = iter_m->current() ;
169 iter_m->next() ;
170}
171
172Vector<uInt> ArrayIndexIteratorAcc::getRows( StorageInitPolicy policy )
173{
174 Block<uInt> v = iter_m->current() ;
175 Int c = isChanged( v ) ;
176// cout << "v=" << Vector<uInt>(IPosition(1,v.nelements()),v.storage(),SHARE) << endl ;
177// cout << "c=" << c << endl ;
178 if ( c > ncol_m-1 ) {
179 pos_m[0] = len_m[ncol_m] ;
180 Int offset = ncol_m - 1 ;
181 while( offset >= 0 && skip_m[offset] )
182 offset-- ;
183 offset++ ;
184// cout << "offset=" << offset << endl ;
185 return Vector<uInt>( pos_m, storage_m.storage()+offset*nrow_m, policy ) ;
186 }
187 Int offset = c - 1 ;
188 while( offset >= 0 && skip_m[offset] )
189 offset-- ;
190 offset++ ;
191// cout << "offset = " << offset << endl ;
192 uInt *base = storage_m.storage() + offset * nrow_m ;
193// cout << "len_m[c+1]=" << len_m[c+1] << endl ;
194// cout << "base=" << Vector<uInt>(IPosition(1,len_m[c+1]),base,SHARE) << endl ;
195 for ( Int i = c ; i < ncol_m ; i++ ) {
196 base = updateStorage( i, base, idxlist_m[i][v[i]] ) ;
197// cout << "len_m[" << i << "]=" << len_m[i] << endl ;
198// cout << "base=" << Vector<uInt>(IPosition(1,len_m[i]),base,SHARE) << endl ;
199 }
200 pos_m[0] = len_m[ncol_m] ;
201// cout << "pos_m=" << pos_m << endl ;
202// cout << "ret=" << Vector<uInt>( pos_m, base, policy ) << endl ;
203 return Vector<uInt>( pos_m, base, policy ) ;
204}
205
206Int ArrayIndexIteratorAcc::isChanged( Block<uInt> &idx )
207{
208 Int i = 0 ;
209 while( i < ncol_m && idx[i] == prev_m[i] )
210 i++ ;
211 return i ;
212}
213
214uInt *ArrayIndexIteratorAcc::updateStorage( Int &icol,
215 uInt *base,
216 uInt &v )
217{
218 // CAUTION:
219 // indexes for storage_m and len_m differ from index for skip_m by 1
220 // (skip_m[0] corresponds to storage_m[1] and len_m[1]) since there
221 // is additional temporary storage at first segment in storage_m
222 uInt *p ;
223 if ( skip_m[icol] ) {
224 // skip update, just update len_m[icol] and pass appropriate pointer
225// cout << "skip " << icol << endl ;
226 p = base ;
227 len_m[icol+1] = len_m[icol] ;
228 }
229 else {
230// cout << "update " << icol << endl ;
231 uInt len = 0 ;
232 p = storage_m.storage() + (icol+1) * nrow_m ;
233 uInt *work = p ;
234 Bool b ;
235 const uInt *arr_p = arr_m.getStorage( b ) ;
236 long offset = 0 ;
237 // warr_p points a first element of (icol)-th column
238 const uInt *warr_p = arr_p + icol * nrow_m ;
239 for ( uInt i = 0 ; i < len_m[icol] ; i++ ) {
240 // increment warr_p by (*(base)-*(base-1))
241 warr_p += *base - offset ;
242 // check if target element is equal to value specified
243 if ( *warr_p == v ) {
244 // then, add current index to storage_m
245 // cout << "add " << *base << endl ;
246 *work = *base ;
247 len++ ;
248 work++ ;
249 }
250 // update offset
251 offset = *base ;
252 // next index
253 base++ ;
254 }
255 arr_m.freeStorage( arr_p, b ) ;
256 len_m[icol+1] = len ;
257 }
258 return p ;
259}
260
261STIdxIter::STIdxIter()
262{
263 iter_m = 0 ;
264}
265
266STIdxIter::STIdxIter( const string &name,
267 const vector<string> &cols )
268{
269}
270
271STIdxIter::STIdxIter( const CountedPtr<Scantable> &s,
272 const vector<string> &cols )
273{
274}
275
276STIdxIter::~STIdxIter()
277{
278 if ( iter_m != 0 )
279 delete iter_m ;
280}
281
282vector<uInt> STIdxIter::tovector( Vector<uInt> v )
283{
284 vector<uInt> ret ;
285 v.tovector( ret ) ;
286 return ret ;
287}
288
289Vector<uInt> STIdxIter::getRows( StorageInitPolicy policy )
290{
291 return iter_m->getRows( policy ) ;
292}
293
294STIdxIterNormal::STIdxIterNormal()
295 : STIdxIter()
296{
297}
298
299STIdxIterNormal::STIdxIterNormal( const string &name,
300 const vector<string> &cols )
301 : STIdxIter( name, cols )
302{
303 Table t( name, Table::Old ) ;
304 init( t, cols ) ;
305}
306
307STIdxIterNormal::STIdxIterNormal( const CountedPtr<Scantable> &s,
308 const vector<string> &cols )
309 : STIdxIter( s, cols )
310{
311 init( s->table(), cols ) ;
312}
313
314STIdxIterNormal::~STIdxIterNormal()
315{
316}
317
318void STIdxIterNormal::init( Table &t,
319 const vector<string> &cols )
320{
321 uInt ncol = cols.size() ;
322 uInt nrow = t.nrow() ;
323 Matrix<uInt> arr( nrow, ncol ) ;
324 ROScalarColumn<uInt> col ;
325 Vector<uInt> v ;
326 for ( uInt i = 0 ; i < ncol ; i++ ) {
327 col.attach( t, cols[i] ) ;
328 v.reference( arr.column( i ) ) ;
329 col.getColumn( v ) ;
330 }
331 iter_m = new ArrayIndexIteratorNormal( arr ) ;
332}
333
334STIdxIterAcc::STIdxIterAcc()
335 : STIdxIter()
336{
337}
338
339STIdxIterAcc::STIdxIterAcc( const string &name,
340 const vector<string> &cols )
341 : STIdxIter( name, cols )
342{
343 Table t( name, Table::Old ) ;
344 init( t, cols ) ;
345}
346
347STIdxIterAcc::STIdxIterAcc( const CountedPtr<Scantable> &s,
348 const vector<string> &cols )
349 : STIdxIter( s, cols )
350{
351 init( s->table(), cols ) ;
352}
353
354STIdxIterAcc::~STIdxIterAcc()
355{
356}
357
358void STIdxIterAcc::init( Table &t,
359 const vector<string> &cols )
360{
361 uInt ncol = cols.size() ;
362 uInt nrow = t.nrow() ;
363 // array shape here is as follows if cols=["BEAMNO","POLNO","IFNO"]:
364 // [[B0,B1,B2,...,BN],
365 // [P0,P1,P2,...,PN],
366 // [I0,I1,I2,...,IN]]
367 // order of internal storage is
368 // [B0,B1,B2,..,BN,P0,P1,P2,...,PN,I0,I1,I2,...,IN]
369 Matrix<uInt> arr( nrow, ncol ) ;
370 Vector<uInt> v ;
371 ROScalarColumn<uInt> col ;
372 for ( uInt i = 0 ; i < ncol ; i++ ) {
373 col.attach( t, cols[i] ) ;
374 v.reference( arr.column( i ) ) ;
375 col.getColumn( v ) ;
376 }
377 iter_m = new ArrayIndexIteratorAcc( arr ) ;
378}
379
380STIdxIterExAcc::STIdxIterExAcc()
381 : STIdxIter(),
382 srctypeid_m( -1 ),
383 srcnameid_m( -1 )
384{
385}
386
387STIdxIterExAcc::STIdxIterExAcc( const string &name,
388 const vector<string> &cols )
389 : STIdxIter( name, cols ),
390 srctypeid_m( -1 ),
391 srcnameid_m( -1 )
392{
393 Table t( name, Table::Old ) ;
394 init( t, cols ) ;
395}
396
397STIdxIterExAcc::STIdxIterExAcc( const CountedPtr<Scantable> &s,
398 const vector<string> &cols )
399 : STIdxIter( s, cols ),
400 srctypeid_m( -1 ),
401 srcnameid_m( -1 )
402{
403 init( s->table(), cols ) ;
404}
405
406STIdxIterExAcc::~STIdxIterExAcc()
407{
408}
409
410void STIdxIterExAcc::init( Table &t,
411 const vector<string> &cols )
412{
413 uInt ncol = cols.size() ;
414 uInt nrow = t.nrow() ;
415 // array shape here is as follows if cols=["BEAMNO","POLNO","IFNO"]:
416 // [[B0,B1,B2,...,BN],
417 // [P0,P1,P2,...,PN],
418 // [I0,I1,I2,...,IN]]
419 // order of internal storage is
420 // [B0,B1,B2,..,BN,P0,P1,P2,...,PN,I0,I1,I2,...,IN]
421 Matrix<uInt> arr( nrow, ncol ) ;
422 Vector<uInt> v ;
423 ROScalarColumn<uInt> col ;
424 ROScalarColumn<String> strCol ;
425 ROScalarColumn<Int> intCol ;
426 for ( uInt i = 0 ; i < ncol ; i++ ) {
427 v.reference( arr.column( i ) ) ;
428 if ( cols[i] == "SRCTYPE" ) {
429 intCol.attach( t, cols[i] ) ;
430 Vector<Int> srctype = intCol.getColumn() ;
431 processIntCol( srctype, v, srctype_m ) ;
432 srctypeid_m = i ;
433 }
434 else if ( cols[i] == "SRCNAME" ) {
435 strCol.attach( t, cols[i] ) ;
436 Vector<String> srcname = strCol.getColumn() ;
437 processStrCol( srcname, v, srcname_m ) ;
438 srcnameid_m = i ;
439 }
440 else {
441 col.attach( t, cols[i] ) ;
442 col.getColumn( v ) ;
443 }
444 }
445 iter_m = new ArrayIndexIteratorAcc( arr ) ;
446}
447
448void STIdxIterExAcc::processIntCol( Vector<Int> &in,
449 Vector<uInt> &out,
450 Block<Int> &val )
451{
452 convertArray( out, in ) ;
453}
454
455void STIdxIterExAcc::processStrCol( Vector<String> &in,
456 Vector<uInt> &out,
457 Block<String> &val )
458{
459 uInt len = in.nelements() ;
460 Vector<String> tmp = in.copy() ;
461 uInt n = genSort( tmp, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates ) ;
462 val.resize( n ) ;
463 for ( uInt i = 0 ; i < n ; i++ ) {
464 val[i] = tmp[i] ;
465// cout << "val[" << i << "]=" << val[i] << endl ;
466 }
467 if ( n == 1 ) {
468 //cout << "n=1" << endl ;
469 out = 0 ;
470 }
471 else if ( n == 2 ) {
472 //cout << "n=2" << endl ;
473 for ( uInt i = 0 ; i < len ; i++ ) {
474 out[i] = (in[i] == val[0]) ? 0 : 1 ;
475 }
476 }
477 else {
478 //cout << "n=" << n << endl ;
479 map<String,uInt> m ;
480 for ( uInt i = 0 ; i < n ; i++ )
481 m[val[i]] = i ;
482 for ( uInt i = 0 ; i < len ; i++ ) {
483 out[i] = m[in[i]] ;
484 }
485 }
486}
487
488Int STIdxIterExAcc::getSrcType()
489{
490 if ( srctypeid_m >= 0 )
491 return (Int)(iter_m->current()[srctypeid_m]) ;
492 else
493 return -999 ;
494}
495
496String STIdxIterExAcc::getSrcName()
497{
498 if ( srcname_m.nelements() > 0 )
499 return srcname_m[iter_m->current()[srcnameid_m]] ;
500 else
501 return "" ;
502}
503
504} // namespace
Note: See TracBrowser for help on using the repository browser.