Commits
Ville Suoranta authored d936bf33dbd Merge
371 371 | Vector<Double> normalizationFactor_p; |
372 372 | VbAvg * vbAvg_p; // [use] |
373 373 | }; |
374 374 | |
375 375 | class VbAvg : public VisBufferImpl2 { |
376 376 | |
377 377 | public: |
378 378 | |
379 379 | friend class MsRowAvg; |
380 380 | |
381 - | VbAvg (const AveragingParameters & averagingParameters, const ViImplementation2 * vi); |
381 + | VbAvg (const AveragingParameters & averagingParameters, ViImplementation2 * vi); |
382 382 | |
383 383 | void accumulate (const VisBuffer2 * vb, const Subchunk & subchunk); |
384 384 | const Cube<Int> & counts () const; |
385 385 | Bool empty () const; |
386 386 | void finalizeBufferFilling (); |
387 387 | void finalizeAverages (); |
388 388 | MsRowAvg * getRow (Int row) const; |
389 389 | MsRowAvg * getRowMutable (Int row); |
390 390 | Bool isComplete () const; |
391 391 | Bool isUsingUvwDistance () const; |
891 891 | void MsRowAvg::setNormalizationFactor(Double normalizationFactor) |
892 892 | { |
893 893 | vbAvg_p->normalizationFactor_p (row ()) = normalizationFactor; |
894 894 | } |
895 895 | |
896 896 | void MsRowAvg::accumulateNormalizationFactor(Double normalizationFactor) |
897 897 | { |
898 898 | vbAvg_p->normalizationFactor_p (row ()) += normalizationFactor; |
899 899 | } |
900 900 | |
901 - | VbAvg::VbAvg (const AveragingParameters & averagingParameters, const ViImplementation2 * vii) |
902 - | : VisBufferImpl2 (VbRekeyable), |
901 + | VbAvg::VbAvg (const AveragingParameters & averagingParameters, ViImplementation2 * vii) |
902 + | : VisBufferImpl2 (vii, VbRekeyable), |
903 903 | averagingInterval_p (averagingParameters.getAveragingInterval ()), |
904 904 | averagingOptions_p (averagingParameters.getOptions()), |
905 905 | averagingVii_p (vii), |
906 906 | bufferToFill_p (0), |
907 907 | complete_p (false), |
908 908 | doing_p (), // all false until determined later on |
909 909 | empty_p (true), |
910 910 | maxTimeDistance_p (averagingParameters.getAveragingInterval() * (0.999)), |
911 911 | // Shrink it just a bit for roundoff |
912 912 | maxUvwDistanceSquared_p (pow(averagingParameters.getMaxUvwDistance(),2)), |
1179 1179 | weightObserved, zeroAccumulation, |
1180 1180 | accumulationParameters.observedOut ()); |
1181 1181 | |
1182 1182 | if (not doing_p.correctedData_p) |
1183 1183 | { |
1184 1184 | // The weight resulting from weighted average is the sum of the weights |
1185 1185 | accumulateElementForCube ( & weightObserved, |
1186 1186 | One, zeroAccumulation, |
1187 1187 | accumulationParameters.weightSpectrumOut ()); |
1188 1188 | } |
1189 - | else |
1190 - | { |
1191 - | // We store the accumulated weight in sigmaSpectrumOut pending of |
1192 - | // - normalization |
1193 - | // - SIGMA = 1/sqrt(WEIGHT) in-place transformation |
1194 - | accumulateElementForCube ( & weightObserved, |
1195 - | One, zeroAccumulation, |
1196 - | accumulationParameters.sigmaSpectrumOut ()); |
1197 - | } |
1189 + | |
1190 + | // This will always create a sigma spectrum column which is not empty. |
1191 + | // This is useful in particular if not doing_p.correctedData_p but doing_p.modelData_p, |
1192 + | // so that modelData can be properly divided by sigmaSpectrumOut in finalizeCubeData |
1193 + | // We store the accumulated weight in sigmaSpectrumOut pending of |
1194 + | // - normalization |
1195 + | // - SIGMA = 1/sqrt(WEIGHT) in-place transformation |
1196 + | accumulateElementForCube (&weightObserved, |
1197 + | One, zeroAccumulation, |
1198 + | accumulationParameters.sigmaSpectrumOut ()); |
1199 + | |
1198 1200 | } |
1199 1201 | |
1200 1202 | // For model data is less clear what to do, what in order to convert to |
1201 1203 | // split we use WEIGHT if averaging CORRECTED_DATA and SIGMA if avg. DATA. |
1202 1204 | // Finally we use WEIGHT by default when averaging MODEL_DATA only |
1203 1205 | if (doing_p.modelData_p) |
1204 1206 | { |
1205 1207 | if (doing_p.correctedData_p) |
1206 1208 | { |
1207 1209 | accumulateElementForCube ( accumulationParameters.modelIn (), |
1475 1477 | struct DividesNonZero : public std::binary_function<L,R,RES> |
1476 1478 | { |
1477 1479 | RES operator() (const L& x, const R& y) const |
1478 1480 | { |
1479 1481 | { return y > 0? RES(x)/y : RES(x); } |
1480 1482 | } |
1481 1483 | }; |
1482 1484 | |
1483 1485 | |
1484 1486 | void |
1485 - | VbAvg::finalizeCubeData (MsRowAvg * msRow) |
1487 + | VbAvg::finalizeCubeData (MsRowAvg * msRowAvg) |
1486 1488 | { |
1487 1489 | // Divide each of the data cubes in use by the sum of the appropriate weights. |
1488 1490 | |
1489 1491 | typedef DividesNonZero <Complex, Float, Complex> DivideOp; |
1490 1492 | DivideOp op; |
1491 1493 | |
1492 1494 | if (doing_p.correctedData_p) |
1493 1495 | { |
1494 - | Matrix<Complex> corrected = msRow->correctedMutable(); |
1495 - | arrayTransformInPlace<Complex, Float, DivideOp > (corrected,msRow->weightSpectrum (), op); |
1496 + | Matrix<Complex> corrected = msRowAvg->correctedMutable(); |
1497 + | arrayTransformInPlace<Complex, Float, DivideOp > (corrected,msRowAvg->weightSpectrum (), op); |
1496 1498 | } |
1497 1499 | |
1498 1500 | if (doing_p.observedData_p) |
1499 1501 | { |
1500 - | Matrix<Complex> observed = msRow->observedMutable(); |
1502 + | Matrix<Complex> observed = msRowAvg->observedMutable(); |
1501 1503 | if (not doing_p.correctedData_p) |
1502 - | { |
1503 - | arrayTransformInPlace<Complex, Float, DivideOp > (observed,msRow->weightSpectrum (), op); |
1504 - | } |
1504 + | arrayTransformInPlace<Complex, Float, DivideOp > (observed,msRowAvg->weightSpectrum (), op); |
1505 1505 | else |
1506 - | { |
1507 - | arrayTransformInPlace<Complex, Float, DivideOp > (observed,msRow->sigmaSpectrum (), op); |
1508 - | } |
1506 + | arrayTransformInPlace<Complex, Float, DivideOp > (observed,msRowAvg->sigmaSpectrum (), op); |
1509 1507 | } |
1510 1508 | |
1511 - | if (doing_p.modelData_p) |
1512 - | { |
1513 - | Matrix<Complex> model = msRow->modelMutable(); |
1509 + | if (doing_p.modelData_p) |
1510 + | { |
1511 + | Matrix<Complex> model = msRowAvg->modelMutable(); |
1514 1512 | |
1515 - | if (doing_p.correctedData_p) |
1516 - | { |
1517 - | arrayTransformInPlace<Complex, Float, DivideOp > (model,msRow->weightSpectrum (), op); |
1518 - | } |
1519 - | else if (doing_p.observedData_p) |
1520 - | { |
1521 - | arrayTransformInPlace<Complex, Float, DivideOp > (model,msRow->sigmaSpectrum (), op); |
1522 - | } |
1523 - | else |
1524 - | { |
1525 - | arrayTransformInPlace<Complex, Int, DivideOp > (model,msRow->counts (), op); |
1526 - | } |
1527 - | } |
1513 + | if (doing_p.correctedData_p) |
1514 + | arrayTransformInPlace<Complex, Float, DivideOp > (model,msRowAvg->weightSpectrum (), op); |
1515 + | else if (doing_p.observedData_p) |
1516 + | arrayTransformInPlace<Complex, Float, DivideOp > (model,msRowAvg->sigmaSpectrum (), op); |
1517 + | else |
1518 + | arrayTransformInPlace<Complex, Int, DivideOp > (model,msRowAvg->counts (), op); |
1519 + | } |
1528 1520 | |
1529 1521 | if (doing_p.floatData_p) |
1530 1522 | { |
1531 1523 | typedef Divides <Float, Float, Float> DivideOpFloat; |
1532 1524 | DivideOpFloat opFloat; |
1533 1525 | |
1534 - | Matrix<Float> visCubeFloat = msRow->singleDishDataMutable(); |
1535 - | arrayTransformInPlace<Float, Float, DivideOpFloat > (visCubeFloat,msRow->weightSpectrum (), opFloat); |
1526 + | Matrix<Float> visCubeFloat = msRowAvg->singleDishDataMutable(); |
1527 + | arrayTransformInPlace<Float, Float, DivideOpFloat > (visCubeFloat,msRowAvg->weightSpectrum (), opFloat); |
1536 1528 | } |
1537 1529 | |
1538 1530 | |
1539 1531 | return; |
1540 1532 | } |
1541 1533 | |
1542 1534 | void |
1543 1535 | VbAvg::finalizeRowData (MsRowAvg * msRow) |
1544 1536 | { |
1545 1537 | Int n = msRow->countsBaseline (); |