Changeset 146
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 -
trunk/src/SDMath.h
r144 r146 45 45 casa::CountedPtr<SDMemTable> quotient(const casa::CountedPtr<SDMemTable>& on, 46 46 const casa::CountedPtr<SDMemTable>& off); 47 47 // 48 48 void multiplyInSitu(SDMemTable* in, casa::Float factor); 49 50 49 casa::CountedPtr<SDMemTable> multiply(const casa::CountedPtr<SDMemTable>& in, 51 52 50 casa::Float factor); 51 // 53 52 casa::CountedPtr<SDMemTable> add(const casa::CountedPtr<SDMemTable>& in, 54 53 casa::Float offset); … … 74 73 enum weightType {NONE,VAR,TSYS}; 75 74 75 void accumulate (casa::Double& timeSum, casa::Double& intSum, casa::Int& nAccum, 76 casa::MaskedArray<casa::Float>& sum, casa::Array<casa::Float>& sumSq, 77 casa::Array<casa::Float>& nPts, casa::Array<casa::Float>& tSysSum, 78 const casa::Array<casa::Float>& tSys, const casa::Array<casa::Float>& nInc, 79 const casa::Vector<casa::Bool>& mask, casa::Double time, casa::Double interval, 80 const casa::Block<casa::CountedPtr<SDMemTable> >& in, 81 casa::uInt iTab, casa::uInt iRow, casa::uInt axis, casa::uInt nAxesSub, 82 casa::Bool useMask, weightType wtType); 83 76 84 void fillSDC (SDContainer& sc, const casa::Array<casa::Bool>& mask, 77 85 const casa::Array<casa::Float>& data, … … 80 88 casa::Double interval, const casa::String& sourceName, 81 89 const casa::Vector<casa::uInt>& freqID); 82 void normalize (casa::MaskedArray<casa::Float>& data,83 const casa::Array<casa::Float>& sumSq,84 const casa::Array<casa::Float>& nPts,85 weightType wtType, casa::Int axis, casa::Int nAxes);86 90 91 SDMemTable* localMultiply (const SDMemTable& in, casa::Float factor); 87 92 93 void normalize (casa::MaskedArray<casa::Float>& data, 94 const casa::Array<casa::Float>& sumSq, 95 const casa::Array<casa::Float>& nPts, 96 weightType wtType, casa::Int axis, casa::Int nAxes); 88 97 }; 89 98 … … 91 100 92 101 #endif 102 103 104 105 106
Note:
See TracChangeset
for help on using the changeset viewer.