Changeset 146 for trunk/src/SDMath.cc
- Timestamp:
- 12/26/04 14:56:06 (19 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/SDMath.cc
r144 r146 268 268 } 269 269 270 // Accumulation step. First get data and deconstruct 271 272 MaskedArray<Float> dataIn(in[iTab]->rowAsMaskedArray(iRow)); 273 Array<Float>& valuesIn = dataIn.getRWArray(); // writable reference 274 const Array<Bool>& maskIn = dataIn.getMask(); // RO reference 275 // 276 if (wtType==NONE) { 277 const MaskedArray<Float> n(nInc,dataIn.getMask()); 278 nPts += n; // Only accumulates where mask==T 279 } else if (wtType==VAR) { 280 281 // We are going to average the data, weighted by the noise for each pol, beam and IF. 282 // So therefore we need to iterate through by spectrum (axis 3) 283 284 VectorIterator<Float> itData(valuesIn, axis); 285 ReadOnlyVectorIterator<Bool> itMask(maskIn, axis); 286 while (!itData.pastEnd()) { 287 288 // Make MaskedArray of Vector, optionally apply OTF mask, and find scaling factor 289 290 if (useMask) { 291 MaskedArray<Float> tmp(itData.vector(),mask&&itMask.vector()); 292 fac = 1.0/variance(tmp); 293 } else { 294 MaskedArray<Float> tmp(itData.vector(),itMask.vector()); 295 fac = 1.0/variance(tmp); 296 } 297 298 // Scale data 299 300 itData.vector() *= fac; // Writes back into 'dataIn' 301 // 302 // Accumulate variance per if/pol/beam averaged over spectrum 303 // This method to get pos2 from itData.pos() is only valid 304 // because the spectral axis is the last one (so we can just 305 // copy the first nAXesSub positions out) 306 307 pos2 = itData.pos().getFirst(nAxesSub); 308 sumSq(pos2) += fac; 309 // 310 itData.next(); 311 itMask.next(); 312 } 313 } else if (wtType==TSYS) { 314 } 315 316 // Accumulate sum of (possibly scaled) data 317 318 sum += dataIn; 319 320 // Accumulate Tsys, time, and interval 321 322 tSysSum += tSys; 323 timeSum += time; 324 intSum += interval; 325 326 // Number of rows in accumulation 327 328 nAccum += 1; 270 // Accumulate 271 272 accumulate (timeSum, intSum, nAccum, sum, sumSq, nPts, tSysSum, 273 tSys, nInc, mask, time, interval, in, iTab, iRow, axis, 274 nAxesSub, useMask, wtType); 275 // 329 276 oldSourceName = sourceName; 330 277 oldFreqID = freqID; … … 433 380 } 434 381 435 void SDMath::multiplyInSitu(SDMemTable* in, Float factor) { 436 SDMemTable* sdmt = new SDMemTable(*in); 437 Table t = sdmt->table(); 438 ArrayColumn<Float> spec(t,"SPECTRA"); 439 for (uInt i=0; i < t.nrow(); i++) { 440 MaskedArray<Float> marr(sdmt->rowAsMaskedArray(i)); 441 marr *= factor; 442 spec.put(i, marr.getArray()); 443 } 444 in = sdmt; 445 delete sdmt;sdmt=0; 446 } 382 383 384 void SDMath::multiplyInSitu(SDMemTable* pIn, Float factor) 385 { 386 SDMemTable* pOut = localMultiply (*pIn, factor); 387 *pIn = *pOut; 388 delete pOut; 389 } 390 447 391 448 392 CountedPtr<SDMemTable> 449 393 SDMath::multiply(const CountedPtr<SDMemTable>& in, Float factor) 450 // 451 // Multiply values by factor 452 // 453 { 454 SDMemTable* sdmt = new SDMemTable(*in); 455 Table t = sdmt->table(); 456 ArrayColumn<Float> spec(t,"SPECTRA"); 457 458 for (uInt i=0; i < t.nrow(); i++) { 459 MaskedArray<Float> marr(sdmt->rowAsMaskedArray(i)); 460 marr *= factor; 461 spec.put(i, marr.getArray()); 462 } 463 return CountedPtr<SDMemTable>(sdmt); 394 { 395 return CountedPtr<SDMemTable>(localMultiply(*in,factor)); 464 396 } 465 397 … … 470 402 // 471 403 { 472 SDMemTable* sdmt = new SDMemTable(*in); 473 404 SDMemTable* sdmt = new SDMemTable(*in,False); 474 405 Table t = sdmt->table(); 475 406 ArrayColumn<Float> spec(t,"SPECTRA"); … … 765 696 } 766 697 698 699 700 // 'private' functions 701 767 702 void SDMath::fillSDC (SDContainer& sc, 768 703 const Array<Bool>& mask, … … 823 758 } 824 759 760 761 void SDMath::accumulate (Double& timeSum, Double& intSum, Int& nAccum, 762 MaskedArray<Float>& sum, Array<Float>& sumSq, 763 Array<Float>& nPts, Array<Float>& tSysSum, 764 const Array<Float>& tSys, const Array<Float>& nInc, 765 const Vector<Bool>& mask, Double time, Double interval, 766 const Block<CountedPtr<SDMemTable> >& in, 767 uInt iTab, uInt iRow, uInt axis, 768 uInt nAxesSub, Bool useMask, 769 weightType wtType) 770 { 771 772 // Get data 773 774 MaskedArray<Float> dataIn(in[iTab]->rowAsMaskedArray(iRow)); 775 Array<Float>& valuesIn = dataIn.getRWArray(); // writable reference 776 const Array<Bool>& maskIn = dataIn.getMask(); // RO reference 777 // 778 if (wtType==NONE) { 779 const MaskedArray<Float> n(nInc,dataIn.getMask()); 780 nPts += n; // Only accumulates where mask==T 781 } else if (wtType==VAR) { 782 783 // We are going to average the data, weighted by the noise for each pol, beam and IF. 784 // So therefore we need to iterate through by spectrum (axis 3) 785 786 VectorIterator<Float> itData(valuesIn, axis); 787 ReadOnlyVectorIterator<Bool> itMask(maskIn, axis); 788 Float fac = 1.0; 789 IPosition pos(nAxesSub,0); 790 // 791 while (!itData.pastEnd()) { 792 793 // Make MaskedArray of Vector, optionally apply OTF mask, and find scaling factor 794 795 if (useMask) { 796 MaskedArray<Float> tmp(itData.vector(),mask&&itMask.vector()); 797 fac = 1.0/variance(tmp); 798 } else { 799 MaskedArray<Float> tmp(itData.vector(),itMask.vector()); 800 fac = 1.0/variance(tmp); 801 } 802 803 // Scale data 804 805 itData.vector() *= fac; // Writes back into 'dataIn' 806 // 807 // Accumulate variance per if/pol/beam averaged over spectrum 808 // This method to get pos2 from itData.pos() is only valid 809 // because the spectral axis is the last one (so we can just 810 // copy the first nAXesSub positions out) 811 812 pos = itData.pos().getFirst(nAxesSub); 813 sumSq(pos) += fac; 814 // 815 itData.next(); 816 itMask.next(); 817 } 818 } else if (wtType==TSYS) { 819 } 820 821 // Accumulate sum of (possibly scaled) data 822 823 sum += dataIn; 824 825 // Accumulate Tsys, time, and interval 826 827 tSysSum += tSys; 828 timeSum += time; 829 intSum += interval; 830 nAccum += 1; 831 } 832 833 SDMemTable* SDMath::localMultiply (const SDMemTable& in, Float factor) 834 { 835 SDMemTable* pOut = new SDMemTable(in,False); 836 const Table& tOut = pOut->table(); 837 ArrayColumn<Float> spec(tOut,"SPECTRA"); 838 // 839 for (uInt i=0; i < tOut.nrow(); i++) { 840 MaskedArray<Float> marr(pOut->rowAsMaskedArray(i)); 841 marr *= factor; 842 spec.put(i, marr.getArray()); 843 } 844 return pOut; 845 } 846 847
Note: See TracChangeset
for help on using the changeset viewer.