Commits

Darrell Schiebel authored 4479ffdd123 Merge
Merge pull request #811 in CASA/casa from feature/CAS-12021 to master

* commit 'ca66d49a6fa02d73bc4e65df260a94ee604add58': Add test for running rerefant on FringeJones solutions Implement re-referencing for FringeJones caltables

code/synthesis/MeasurementComponents/FringeJones.cc

Modified
22 22 //# National Radio Astronomy Observatory
23 23 //# 520 Edgemont Road
24 24 //# Charlottesville, VA 22903-2475 USA
25 25 //#
26 26
27 27 #include <synthesis/MeasurementComponents/FringeJones.h>
28 28
29 29 #include <msvis/MSVis/VisBuffer.h>
30 30 #include <msvis/MSVis/VisBuffAccumulator.h>
31 31 #include <ms/MeasurementSets/MSColumns.h>
32 +#include <synthesis/CalTables/CTIter.h>
32 33 #include <synthesis/MeasurementEquations/VisEquation.h> // *
33 34 #include <synthesis/MeasurementComponents/SolveDataBuffer.h>
35 +#include <synthesis/MeasurementComponents/MSMetaInfoForCal.h>
34 36 #include <lattices/Lattices/ArrayLattice.h>
35 37 #include <lattices/LatticeMath/LatticeFFT.h>
36 38 #include <scimath/Mathematics/FFTServer.h>
37 39
38 40 #include <casa/Arrays/ArrayMath.h>
39 41 #include <casa/Arrays/MatrixMath.h>
40 42 #include <casa/Arrays/ArrayLogical.h>
41 43 #include <casa/BasicSL/String.h>
42 44 #include <casa/Utilities/Assert.h>
43 45 #include <casa/Exceptions/Error.h>
2294 2296 }
2295 2297 }
2296 2298 }
2297 2299
2298 2300 void
2299 2301 FringeJones::solveOneVB(const VisBuffer&) {
2300 2302 throw(AipsError("VisBuffer interface not supported!"));
2301 2303 }
2302 2304
2303 2305
2306 +void FringeJones::globalPostSolveTinker() {
2307 +
2308 + // Re-reference the phase, if requested
2309 + if (refantlist()(0)>-1) applyRefAnt();
2310 +}
2311 +
2312 +void FringeJones::applyRefAnt() {
2313 +
2314 + // TBD:
2315 + // 1. Synchronize refant changes on par axis
2316 + // 2. Implement minimum mean deviation algorithm
2317 +
2318 + if (refantlist()(0)<0)
2319 + throw(AipsError("No refant specified."));
2320 +
2321 + Int nUserRefant=refantlist().nelements();
2322 +
2323 + // Get the preferred refant names from the MS
2324 + String refantName(msmc().antennaName(refantlist()(0)));
2325 + if (nUserRefant>1) {
2326 + refantName+=" (";
2327 + for (Int i=1;i<nUserRefant;++i) {
2328 + refantName+=msmc().antennaName(refantlist()(i));
2329 + if (i<nUserRefant-1) refantName+=",";
2330 + }
2331 + refantName+=")";
2332 + }
2333 +
2334 + logSink() << "Applying refant: " << refantName
2335 + << " refantmode = " << refantmode();
2336 + if (refantmode()=="flex")
2337 + logSink() << " (hold alternate refants' phase constant) when refant flagged";
2338 + if (refantmode()=="strict")
2339 + logSink() << " (flag all antennas when refant flagged)";
2340 + logSink() << LogIO::POST;
2341 +
2342 + // Generate a prioritized refant choice list
2343 + // The first entry in this list is the user's primary refant,
2344 + // the second entry is the refant used on the previous interval,
2345 + // and the rest is a prioritized list of alternate refants,
2346 + // starting with the user's secondary (if provided) refants,
2347 + // followed by the rest of the array, in distance order. This
2348 + // makes the priorities correct all the time, and prevents
2349 + // a semi-stochastic alternation (by preferring the last-used
2350 + // alternate, even if nominally higher-priority refants become
2351 + // available)
2352 +
2353 +
2354 + // Extract antenna positions
2355 + Matrix<Double> xyz;
2356 + if (msName()!="<noms>") {
2357 + MeasurementSet ms(msName());
2358 + ROMSAntennaColumns msant(ms.antenna());
2359 + msant.position().getColumn(xyz);
2360 + }
2361 + else {
2362 + // TBD RO*
2363 + CTColumns ctcol(*ct_);
2364 + CTAntennaColumns& antcol(ctcol.antenna());
2365 + antcol.position().getColumn(xyz);
2366 + }
2367 +
2368 + // Calculate (squared) antenna distances, relative
2369 + // to last preferred antenna
2370 + Vector<Double> dist2(xyz.ncolumn(),0.0);
2371 + for (Int i=0;i<3;++i) {
2372 + Vector<Double> row=xyz.row(i);
2373 + row-=row(refantlist()(nUserRefant-1));
2374 + dist2+=square(row);
2375 + }
2376 + // Move preferred antennas to a large distance
2377 + for (Int i=0;i<nUserRefant;++i)
2378 + dist2(refantlist()(i))=DBL_MAX;
2379 +
2380 + // Generated sorted index
2381 + Vector<uInt> ord;
2382 + genSort(ord,dist2);
2383 +
2384 + // Assemble the whole choices list
2385 + Int nchoices=nUserRefant+1+ord.nelements();
2386 + Vector<Int> refantchoices(nchoices,0);
2387 + Vector<Int> r(refantchoices(IPosition(1,nUserRefant+1),IPosition(1,refantchoices.nelements()-1)));
2388 + convertArray(r,ord);
2389 +
2390 + // set first two to primary preferred refant
2391 + refantchoices(0)=refantchoices(1)=refantlist()(0);
2392 +
2393 + // set user's secondary refants (if any)
2394 + if (nUserRefant>1)
2395 + refantchoices(IPosition(1,2),IPosition(1,nUserRefant))=
2396 + refantlist()(IPosition(1,1),IPosition(1,nUserRefant-1));
2397 +
2398 + //cout << "refantchoices = " << refantchoices << endl;
2399 +
2400 +
2401 + if (refantmode()=="strict") {
2402 + nchoices=1;
2403 + refantchoices.resize(1,True);
2404 + }
2405 +
2406 + Vector<Int> nPol(nSpw(),nPar()); // TBD:or 1, if data was single pol
2407 +
2408 + if (nPar()==6) {
2409 + // Verify that 2nd poln has unflagged solutions, PER SPW
2410 + ROCTMainColumns ctmc(*ct_);
2411 +
2412 + Block<String> cols(1);
2413 + cols[0]="SPECTRAL_WINDOW_ID";
2414 + CTIter ctiter(*ct_,cols);
2415 + Cube<Bool> fl;
2416 +
2417 + while (!ctiter.pastEnd()) {
2418 +
2419 + Int ispw=ctiter.thisSpw();
2420 + fl.assign(ctiter.flag());
2421 +
2422 + IPosition blc(3,0,0,0), trc(fl.shape());
2423 + trc-=1; trc(0)=blc(0)=1;
2424 +
2425 + // cout << "ispw = " << ispw << " nfalse(fl(1,:,:)) = " << nfalse(fl(blc,trc)) << endl;
2426 +
2427 + // If there are no unflagged solutions in 2nd pol,
2428 + // avoid it in refant calculations
2429 + if (nfalse(fl(blc,trc))==0)
2430 + nPol(ispw)=1;
2431 +
2432 + ctiter.next();
2433 + }
2434 + }
2435 + // cout << "nPol = " << nPol << endl;
2436 +
2437 + Bool usedaltrefant(false);
2438 + Int currrefant(refantchoices(0)), lastrefant(-1);
2439 +
2440 + Block<String> cols(2);
2441 + cols[0]="SPECTRAL_WINDOW_ID";
2442 + cols[1]="TIME";
2443 + CTIter ctiter(*ct_,cols);
2444 +
2445 + // Arrays to hold per-timestamp solutions
2446 + Cube<Float> solA, solB;
2447 + Cube<Bool> flA, flB;
2448 + Vector<Int> ant1A, ant1B, ant2B;
2449 + Matrix<Complex> refPhsr; // the reference phasor [npol,nchan]
2450 + Int lastspw(-1);
2451 + Bool first(true);
2452 + while (!ctiter.pastEnd()) {
2453 + Int ispw=ctiter.thisSpw();
2454 + if (ispw!=lastspw) first=true; // spw changed, start over
2455 +
2456 + // Read in the current sol, fl, ant1:
2457 + solB.assign(ctiter.fparam());
2458 + flB.assign(ctiter.flag());
2459 + ant1B.assign(ctiter.antenna1());
2460 + ant2B.assign(ctiter.antenna2());
2461 +
2462 + // First time thru, 'previous' solution same as 'current'
2463 + if (first) {
2464 + solA.reference(solB);
2465 + flA.reference(flB);
2466 + ant1A.reference(ant1B);
2467 + }
2468 + IPosition shB(solB.shape());
2469 + IPosition shA(solA.shape());
2470 +
2471 + // Find a good refant at this time
2472 + // A good refant is one that is unflagged in all polarizations
2473 + // in the current(B) and previous(A) intervals (so they can be connected)
2474 + Int irefA(0),irefB(0); // index on 3rd axis of solution arrays
2475 + Int ichoice(0); // index in refantchoicelist
2476 + Bool found(false);
2477 + IPosition blcA(3,0,0,0),trcA(shA),blcB(3,0,0,0),trcB(shB);
2478 + trcA-=1; trcA(0)=trcA(2)=0;
2479 + trcB-=1; trcB(0)=trcB(2)=0;
2480 + ichoice=0;
2481 + while (!found && ichoice<nchoices) {
2482 + // Find index of current refant choice
2483 + irefA=irefB=0;
2484 + while (ant1A(irefA)!=refantchoices(ichoice) && irefA<shA(2)) ++irefA;
2485 + while (ant1B(irefB)!=refantchoices(ichoice) && irefB<shB(2)) ++irefB;
2486 +
2487 + if (irefA<shA(2) && irefB<shB(2)) {
2488 +
2489 + // cout << " Trial irefA,irefB: " << irefA << "," << irefB
2490 + // << " Ants=" << ant1A(irefA) << "," << ant1B(irefB) << endl;
2491 +
2492 + blcA(2)=trcA(2)=irefA;
2493 + blcB(2)=trcB(2)=irefB;
2494 + found=true; // maybe
2495 + for (Int ipol=0;ipol<nPol(ispw);++ipol) {
2496 + blcA(0)=trcA(0)=blcB(0)=trcB(0)=ipol;
2497 + found &= (nfalse(flA(blcA,trcA))>0); // previous interval
2498 + found &= (nfalse(flB(blcB,trcB))>0); // current interval
2499 + }
2500 + }
2501 + else
2502 + // irefA or irefB out-of-range
2503 + found=false; // Just to be sure
2504 +
2505 + if (!found) ++ichoice; // try next choice next round
2506 +
2507 + }
2508 +
2509 + if (found) {
2510 + // at this point, irefA/irefB point to a good refant
2511 +
2512 + // Keep track
2513 + usedaltrefant|=(ichoice>0);
2514 + currrefant=refantchoices(ichoice);
2515 + refantchoices(1)=currrefant; // 2nd priorty next time
2516 +
2517 + // cout << " currrefant = " << currrefant << " (" << ichoice << ")" << endl;
2518 +
2519 + // cout << " Final irefA,irefB: " << irefA << "," << irefB
2520 + // << " Ants=" << ant1A(irefA) << "," << ant1B(irefB) << endl;
2521 +
2522 +
2523 + // Only report if using an alternate refant
2524 + if (currrefant!=lastrefant && ichoice>0) {
2525 + logSink()
2526 + << "At "
2527 + << MVTime(ctiter.thisTime()/C::day).string(MVTime::YMD,7)
2528 + << " ("
2529 + << "Spw=" << ctiter.thisSpw()
2530 + << ", Fld=" << ctiter.thisField()
2531 + << ")"
2532 + << ", using refant " << msmc().antennaName(currrefant)
2533 + << " (id=" << currrefant
2534 + << ")" << " (alternate)"
2535 + << LogIO::POST;
2536 + }
2537 +
2538 + // Form reference phasor [nPar,nChan]
2539 + Matrix<Float> rA,rB;
2540 + Matrix<Bool> rflA,rflB;
2541 + rB.assign(solB.xyPlane(irefB));
2542 + rflB.assign(flB.xyPlane(irefB));
2543 +
2544 + if (!first) {
2545 + // Get and condition previous phasor for the current refant
2546 + rA.assign(solA.xyPlane(irefA));
2547 + rflA.assign(flA.xyPlane(irefA));
2548 + rB-=rA;
2549 +
2550 + // Accumulate flags
2551 + rflB&=rflA;
2552 + }
2553 +
2554 + // cout << " rB = " << rB << endl;
2555 + // cout << boolalpha << " rflB = " << rflB << endl;
2556 + // TBD: fillChanGaps?
2557 +
2558 + // Now apply reference phasor to all antennas
2559 + Matrix<Float> thissol;
2560 + for (Int iant=0;iant<shB(2);++iant) {
2561 + thissol.reference(solB.xyPlane(iant));
2562 + thissol-=rB;
2563 + }
2564 +
2565 + // Set refant, so we can put it back
2566 + ant2B=currrefant;
2567 +
2568 + // put back referenced solutions
2569 + ctiter.setfparam(solB);
2570 + ctiter.setantenna2(ant2B);
2571 +
2572 + // Next time thru, solB is previous
2573 + solA.reference(solB);
2574 + flA.reference(flB);
2575 + ant1A.reference(ant1B);
2576 + solB.resize(); // (break references)
2577 + flB.resize();
2578 + ant1B.resize();
2579 +
2580 + lastrefant=currrefant;
2581 + first=false; // avoid first-pass stuff from now on
2582 +
2583 + } // found
2584 + else {
2585 + logSink()
2586 + << "At "
2587 + << MVTime(ctiter.thisTime()/C::day).string(MVTime::YMD,7)
2588 + << " ("
2589 + << "Spw=" << ctiter.thisSpw()
2590 + << ", Fld=" << ctiter.thisField()
2591 + << ")"
2592 + << ", refant (id=" << currrefant
2593 + << ") was flagged; flagging all antennas strictly."
2594 + << LogIO::POST;
2595 + // Flag all solutions in this interval
2596 + flB.set(True);
2597 + ctiter.setflag(flB);
2598 + }
2599 +
2600 + // advance to the next interval
2601 + lastspw=ispw;
2602 + ctiter.next();
2603 + }
2604 +
2605 + if (usedaltrefant)
2606 + logSink() << LogIO::NORMAL
2607 + << " NB: An alternate refant was used at least once to maintain" << endl
2608 + << " phase continuity where the user's preferred refant drops out." << endl
2609 + << " Alternate refants are held constant in phase (_not_ zeroed)" << endl
2610 + << " during these periods, and the preferred refant may return at" << endl
2611 + << " a non-zero phase. This is generally harmless."
2612 + << LogIO::POST;
2613 +
2614 + return;
2615 +
2616 +}
2617 +
2304 2618 } //# NAMESPACE CASA - END
2305 2619
2306 2620
2307 2621 /*
2308 2622 Example of use in at Python level:
2309 2623
2310 2624 fringefit(vis="n14c2.ms", caltable="fail.fj", field="",spw="1",intent="",
2311 2625 selectdata=True, timerange="", antenna="", scan="5", observation="",
2312 2626 msselect="", solint="inf", refant="EF", minsnr=3.0, append=False,
2313 2627 gaintable=['n14c2.gcal'], parang=False)

Everything looks good. We'll let you know here if there's anything you should know about.

Add shortcut