Commits

George Moellenbrock authored 76c9b92fdbe
Introduce L1 and R solving capabilities (gaincal only).
No tags

code/synthesis/MeasurementComponents/test/tGJones_GT.cc

Modified
30 30 #include <casa/iostream.h>
31 31 #include <msvis/MSVis/VisBuffer2.h>
32 32 #include <msvis/MSVis/SimpleSimVi2.h>
33 33 //#include <synthesis/MeasurementComponents/SolveDataBuffer.h>
34 34 #include <casa/Arrays/ArrayMath.h>
35 35 #include <casa/OS/Timer.h>
36 36 #include <synthesis/MeasurementComponents/StandardVisCal.h>
37 37 #include <synthesis/MeasurementComponents/DJones.h>
38 38 #include <synthesis/MeasurementComponents/KJones.h>
39 39 #include <synthesis/MeasurementComponents/MSMetaInfoForCal.h>
40 +#include <synthesis/MeasurementEquations/VisEquation.h>
41 +#include <synthesis/MeasurementComponents/VisCalSolver2.h>
40 42
41 43 #include <gtest/gtest.h>
42 44
45 +#include "VisCalTestBase_GT.h"
46 +
43 47 #define SHOWSTATE false
44 48 using namespace casacore;
45 49 using namespace casa;
46 50 using namespace casa::vi;
47 51
48 52 class VisCalTest : public ::testing::Test {
49 53
50 54 public:
51 55
52 56 VisCalTest() :
345 349 //cout << "dblph=" << blph << endl;
346 350
347 351 ASSERT_TRUE(allNearAbs(blph,0.0f,5e-5));
348 352
349 353 delete G;
350 354
351 355 Table::deleteTable("refant_flex.G");
352 356 Table::deleteTable("refant_strict.G");
353 357
354 358 }
359 +
360 +
361 +
362 +class GJonesSolveTest : public VisCalTestBase {
363 +
364 +public:
365 +
366 + Cube<Complex> g;
367 +
368 + GJonesSolveTest() :
369 + VisCalTestBase(1,1,1,10,4,1,1,false,"circ"),
370 + g(2,1,nAnt,Complex(1.0))
371 + {
372 + /*
373 + for (Int iant=0;iant<nAnt;++iant) {
374 + Float P(iant*C::pi/2);
375 + //g(0,0,iant)=(1.0+iant*0.01)*Complex(cos(P),sin(P));
376 + //g(1,0,iant)=(1.0-iant*0.01)*Complex(cos(P),sin(-P));
377 + g(0,0,iant)=Complex(cos(P),sin(P));
378 + g(1,0,iant)=Complex(cos(P),sin(-P));
379 + }
380 + */
381 +
382 + summary("GJonesSolveTest");
383 + cout << "g = " << g << endl;
384 + }
385 +
386 +};
387 +
388 +
389 +TEST_F(GJonesSolveTest, Test1) {
390 +
391 + // Apply-able Xf
392 + GJones Gapp(msmc);
393 + Gapp.setApply();
394 +
395 + for (Int ispw=0;ispw<nSpw;++ispw) {
396 + Gapp.setMeta(0,1,0.0,
397 + ispw,ssvp.freqs(ispw),
398 + nChan);
399 + Gapp.sizeApplyParCurrSpw(nChan);
400 +
401 + Gapp.setApplyParCurrSpw(g,true,false); // correct below will corrupt
402 + }
403 + //Gapp.state();
404 +
405 +
406 + Record solvePar;
407 + solvePar.define("table",String("test.G"));
408 + solvePar.define("solint",String("inf"));
409 + solvePar.define("combine",String(""));
410 + solvePar.define("minblperant",3);
411 + solvePar.define("solmode",String(""));
412 + //Vector<Double> rmsth(std::vector<Float>({3.0,2.5}));
413 + //solvePar.define("rmsthresh",rmsth);
414 + Vector<Int> refant(1,0); solvePar.define("refant",refant);
415 +
416 + cout << "solvePar= " << solvePar << endl;
417 +
418 + Vector<String> solmodes(4,"");
419 + solmodes(0)="";
420 + solmodes(1)="L1";
421 + solmodes(2)="R";
422 + solmodes(3)="L1R";
423 +
424 +
425 + cout << "-----------" << endl;
426 +
427 + for (Int ism=0;ism<4;++ism) {
428 +
429 + cout << "solmode=" << solmodes(ism) << endl;
430 +
431 + solvePar.define("solmode",solmodes(ism));
432 +
433 + GJones Gsol(msmc);
434 + Gsol.setSolve(solvePar);
435 + // Gsol.state();
436 + Gsol.createMemCalTable2();
437 +
438 +
439 + for (vi2.originChunks();vi2.moreChunks();vi2.nextChunk()) {
440 + for (vi2.origin();vi2.more();vi2.next()) {
441 +
442 + Int ispw=vb2->spectralWindows()(0);
443 + Int obsid(vb2->observationId()(0));
444 + Int scan(vb2->scan()(0));
445 + Double timestamp(vb2->time()(0));
446 + Int fldid(vb2->fieldId()(0));
447 + Vector<Double> freqs(vb2->getFrequencies(0));
448 + Vector<Int> a1(vb2->antenna1());
449 + Vector<Int> a2(vb2->antenna2());
450 +
451 + vb2->resetWeightsUsingSigma();
452 +
453 + Cube<Complex> vC(vb2->visCube());
454 + Float onedeg(C::pi/180.0);
455 + vC*=Complex(cos(onedeg),sin(onedeg));
456 +
457 + //vC(0,0,6)*=Complex(0.0,sin(C::pi/2.0));
458 + //vC(0,0,26)*=Complex(-1.0,0.0);
459 + vC(0,0,6)*=0.5;
460 + vC(0,0,11)*=2.0;
461 + vC(0,0,17)*=2.0;
462 + vC(0,0,18)*=0.5;
463 + vC(0,0,35)*=1.9;
464 + vC(0,0,44)*=1.05;
465 +
466 +
467 + vb2->setVisCubeCorrected(vC);
468 + vb2->setFlagCube(vb2->flagCube());
469 +
470 + Gapp.setMeta(obsid,scan,timestamp,
471 + ispw,freqs,
472 + fldid);
473 + Gapp.correct2(*vb2,false,false,false); // (trial?,doWtSp?,dosync?)
474 +
475 + /*
476 + Int nrow=vb2->nRows();
477 + Double t0(86400.0*floor(timestamp/86400.0));
478 + for (Int irow=0;irow<nrow;++irow) {
479 + cout << iros << " time=" << timestamp-t0
480 + << " BL=" << a1[irow] << "-" << a2[irow]
481 + << " vCC=" << vb2->visCubeCorrected()(Slice(),Slice(0),Slice(irow)).nonDegenerate() << endl;
482 + }
483 + */
484 + SDBList sdbs;
485 + sdbs.add(*vb2);
486 + sdbs.enforceSolveWeights(Gsol.phandonly());
487 +
488 + // Use syncSolveMeta here??
489 +
490 + // Setup meta & sizes for the solve
491 + Gsol.setMeta(sdbs.aggregateObsId(),
492 + sdbs.aggregateScan(),
493 + sdbs.aggregateTime(),
494 + sdbs.aggregateSpw(),
495 + sdbs.freqs(),
496 + sdbs.aggregateFld());
497 +
498 + Gsol.setOrVerifyCTFrequencies(sdbs.aggregateSpw());
499 + Gsol.sizeSolveParCurrSpw(sdbs.nChannels());
500 + //Gsol.state();
501 +
502 +
503 + Gsol.guessPar(sdbs);
504 + //cout << "solvePar() = " << Gsol.solveCPar() << endl;
505 +
506 + // Munge LL flags _after_ guess, so L solutions are "correct"
507 + if (solmodes(ism)!="") {
508 + Cube<Bool> flc(sdbs(0).flagCube());
509 + flc(3,0,6)=true;
510 + flc(3,0,11)=true;
511 + flc(3,0,17)=true;
512 + flc(3,0,18)=true;
513 + flc(3,0,35)=true;
514 + flc(3,0,44)=true;
515 + }
516 +
517 + VisEquation ve;
518 + ve.setsolve(Gsol);
519 + VisCalSolver2 vcs(Gsol.solmode(),Gsol.rmsthresh());
520 +
521 + Bool good=vcs.solve(ve,Gsol,sdbs);
522 + //cout << "good=" << good << endl;
523 + ASSERT_TRUE(good);
524 +
525 + // Gsol.applyRefAnt();
526 +
527 + //cout << "solvePar() = " << Gsol.solveCPar() << endl;
528 + //cout << "solveParErr() = " << Gsol.solveParErr() << endl;
529 +
530 + // Form phase and amp
531 + Cube<Complex> soln;
532 + soln.assign(Gsol.solveCPar());
533 + for (Int i=0;i<2;++i) {
534 + Complex Cph=soln(i,0,0);
535 + Cph/=abs(Cph);
536 + Cube<Complex> solnsl(soln(Slice(i,1,1),Slice(),Slice()));
537 + solnsl/=Cph;
538 + }
539 + Matrix<Float> Gamp(amplitude(soln).nonDegenerate());
540 + Matrix<Float> Gpha(phase(soln).nonDegenerate()*180.0/C::pi);
541 +
542 + //cout << "amplitudes=" << amplitude(Gsol.solveCPar()) << endl;
543 + //cout << "phases =" << phase(soln)*180.0/C::pi << endl;
544 +
545 + //cout << "amplitudes=" << Gamp << endl;
546 + //cout << "phases =" << Gpha << endl;
547 +
548 + //cout << "Residuals: " << amplitude(sdbs(0).residuals()) << endl;
549 +
550 + // Form R/L ratio for testing
551 + Vector<Complex> L(soln(Slice(1),Slice(),Slice()).nonDegenerate());
552 + Vector<Complex> RoL(soln(Slice(0),Slice(),Slice()).nonDegenerate());
553 + RoL/=L;
554 + RoL-=Complex(1.0);
555 + Vector<Float> Arol(amplitude(RoL));
556 + Float maxArol(max(Arol));
557 +
558 + //cout << " R/L = " << Arol
559 + //<< " max = " << maxArol
560 + // << endl;
561 +
562 +
563 + if (solmodes(ism)=="") {
564 + // This is the traditional solution
565 +
566 + // Expecting L solutions to be correct, R solutions to be no worse than 10% off
567 +
568 + // Find peak absolute amp resid
569 + Gamp-=1.0f;
570 + Gamp=abs(Gamp);
571 + Vector<Float> maxAresid(partialMaxs(Gamp,IPosition(1,1)));
572 +
573 + //cout << "Gamp_resid = " << Gamp << endl;
574 + cout << "max(Gamp_resid) = " << maxAresid << endl;
575 + ASSERT_TRUE(maxAresid(0)<0.1);
576 + ASSERT_TRUE(maxAresid(1)<1e-4);
577 +
578 + // Find peak absolute phase resid (deg)
579 + for (Int iant=1;iant<nAnt;++iant) {
580 + Vector<Float> Piant(Gpha(Slice(),Slice(iant,1,1)));
581 + Piant+=Float(iant*0.2);
582 + }
583 + Vector<Float> maxPresid(partialMaxs(Gpha,IPosition(1,1)));
584 +
585 + //cout << "Gpha_resid = " << Gpha << endl;
586 + cout << "max(Gpha_resid) = " << maxPresid << endl;
587 + ASSERT_TRUE(maxPresid(0)<0.02);
588 + ASSERT_TRUE(maxPresid(1)<1e-5);
589 +
590 + // Test max R/L
591 + cout << "maxArol = " << maxArol << endl;
592 + ASSERT_TRUE(maxArol<0.1f);
593 +
594 +
595 + }
596 + if (solmodes(ism).contains("R")) {
597 + Matrix<Bool> wFl(sdbs(0).const_workingFlagCube().nonDegenerate());
598 +
599 + //cout << "flags=" << boolalpha << wFl << endl;
600 + ASSERT_TRUE(ntrue(wFl)==12);
601 +
602 + ASSERT_TRUE(wFl(0,6));
603 + ASSERT_TRUE(wFl(0,11));
604 + ASSERT_TRUE(wFl(0,17));
605 + ASSERT_TRUE(wFl(0,18));
606 + ASSERT_TRUE(wFl(0,35));
607 + ASSERT_TRUE(wFl(0,44));
608 +
609 + if (solmodes(ism).contains("L1")) {
610 + Matrix<Float> wWt(sdbs(0).const_workingWtSpec().nonDegenerate());
611 +
612 + //cout << "wts=" << wWt << endl;
613 + ASSERT_EQ(wWt(0,6),0.0f);
614 + ASSERT_EQ(wWt(0,11),0.0f);
615 + ASSERT_EQ(wWt(0,17),0.0f);
616 + ASSERT_EQ(wWt(0,18),0.0f);
617 + ASSERT_EQ(wWt(0,35),0.0f);
618 + ASSERT_EQ(wWt(0,44),0.0f);
619 + }
620 +
621 + }
622 +
623 + // Test max R/L by mode
624 + if (solmodes(ism)=="L1") {
625 + cout << "maxArol = " << maxArol << endl;
626 + ASSERT_TRUE(maxArol<3e-3);
627 + }
628 + if (solmodes(ism)=="R") {
629 + cout << "maxArol = " << maxArol << endl;
630 + ASSERT_TRUE(maxArol<1e-6);
631 + }
632 + if (solmodes(ism)=="L1R") {
633 + cout << "maxArol = " << maxArol << endl;
634 + ASSERT_TRUE(maxArol<2e-7);
635 + }
636 +
637 + }
638 + }
639 +
640 + cout << "-----------" << endl;
641 +
642 + } // solmodes
643 +
644 +
645 +}

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

Add shortcut