Commits
Pam Ford authored 5a7613aeeee Merge
39 39 | // Test program for FringeJones |
40 40 | // </summary> |
41 41 | |
42 42 | |
43 43 | using namespace std; |
44 44 | using namespace casa; |
45 45 | using namespace casacore; |
46 46 | using namespace casa::vi; |
47 47 | |
48 48 | |
49 - | |
49 + | |
50 50 | |
51 51 | int main(int argc, char **argv) { |
52 52 | ::testing::InitGoogleTest(&argc, argv); |
53 53 | return RUN_ALL_TESTS(); |
54 54 | } |
55 55 | |
56 56 | class DelayRateFFTTest : public ::testing::Test { |
57 57 | |
58 58 | public: |
59 59 | virtual Array<Complex> appdel(Int nchan, Int nt, Float f0, Float df, Float dt, Float delay, Float fringerate) { |
78 78 | Float rate(4.0e-13); |
79 79 | SDBList s = SDBList(); |
80 80 | |
81 81 | Array<Complex> Vobs0(IPosition(4, 1, nelem, nt, nchan)); |
82 82 | const Array<Complex>& V2 = Vobs0(Slicer(IPosition(4, 0, nelem-1, 0, 0), |
83 83 | IPosition(4, 0, nelem-1, nt-1, nchan-1), |
84 84 | IPosition(4, 1, 1, 1, 1), |
85 85 | Slicer::endIsLast)).nonDegenerate(); |
86 86 | const Array<Complex>& V3 = this->appdel(nchan, nt, f0, df, dt, delay, rate); |
87 87 | V2.nonDegenerate() = V3; |
88 + | |
89 + | Array<Double> delayWindow(IPosition(1, 2)); |
90 + | Array<Double> rateWindow(IPosition(1, 2)); |
91 + | delayWindow(IPosition(1, 0)) = -100.0; |
92 + | delayWindow(IPosition(1, 1)) = +100.0; |
93 + | rateWindow(IPosition(1, 0)) = -100.0; |
94 + | rateWindow(IPosition(1, 1)) = +100.0; |
95 + | |
96 + | IPosition ds = delayWindow.shape(); |
88 97 | |
89 - | DelayRateFFT drfft0(Vobs0, nPadfactor, f0, df, dt, s); |
98 + | DelayRateFFT drfft0(Vobs0, nPadfactor, f0, df, dt, s, delayWindow, rateWindow); |
90 99 | drfft0.FFT(); |
91 100 | drfft0.searchPeak(); |
92 101 | Float rate_resn = 1.0f/(nPadfactor*nt*dt*1e9*f0); |
93 102 | Float delay_resn = 1.0f/(nPadfactor*nchan*df); |
94 103 | |
95 104 | if (FRINGEJONES_TEST_VERBOSE) { |
96 105 | cerr << boolalpha; |
97 106 | cerr << "param = " << drfft0.param() << endl; |
98 107 | cerr << "delay = " << drfft0.delay()(0,1) << endl; |
99 - | cerr << "delay resolution =" << delay_resn << endl; |
100 - | cerr << "scaled errror in delay = " << drfft0.delay()(0,1) - delay << endl; |
108 + | cerr << "delay resolution = " << delay_resn << endl; |
109 + | cerr << "scaled error in delay = " << (drfft0.delay()(0, 1) - delay)/delay_resn << endl; |
101 110 | cerr << "rate = " << drfft0.rate()(0, 1) << endl; |
102 - | cerr << "rate resolution =" << rate_resn << endl; |
111 + | cerr << "rate resolution = " << rate_resn << endl; |
103 112 | cerr << "scaled error in rate = " << (drfft0.rate()(0, 1) - rate)/rate_resn << endl; |
104 113 | } |
105 114 | ASSERT_TRUE(allNearAbs(drfft0.delay()(0,1), delay, 0.5*delay_resn)); |
106 115 | ASSERT_TRUE(allNearAbs(drfft0.rate()(0,1), rate, 0.5*rate_resn)); |
107 116 | } |
108 117 | |
109 118 | |
119 + | TEST_F(DelayRateFFTTest, BasicDelayRateFFTTest2) { |
120 + | Int nPadfactor(4); |
121 + | Int nchan(32), nt(20), nelem(2); |
122 + | Float df(0.1), dt(2.0);; |
123 + | Float f0(9.0); |
124 + | Float delay(-2.8); |
125 + | Float rate(4.0e-13); |
126 + | SDBList s = SDBList(); |
127 + | |
128 + | Array<Complex> Vobs0(IPosition(4, 1, nelem, nt, nchan)); |
129 + | const Array<Complex>& V2 = Vobs0(Slicer(IPosition(4, 0, nelem-1, 0, 0), |
130 + | IPosition(4, 0, nelem-1, nt-1, nchan-1), |
131 + | IPosition(4, 1, 1, 1, 1), |
132 + | Slicer::endIsLast)).nonDegenerate(); |
133 + | const Array<Complex>& V3 = this->appdel(nchan, nt, f0, df, dt, delay, rate); |
134 + | V2.nonDegenerate() = V3; |
135 + | |
136 + | Array<Double> delayWindow(IPosition(1, 2)); |
137 + | Array<Double> rateWindow(IPosition(1, 2)); |
138 + | Float minDelay = -1.0; |
139 + | Float maxDelay = +1.0; |
140 + | delayWindow(IPosition(1, 0)) = minDelay; |
141 + | delayWindow(IPosition(1, 1)) = maxDelay; |
142 + | rateWindow(IPosition(1, 0)) = -100.0; |
143 + | rateWindow(IPosition(1, 1)) = +100.0; |
144 + | |
145 + | IPosition ds = delayWindow.shape(); |
146 + | |
147 + | DelayRateFFT drfft0(Vobs0, nPadfactor, f0, df, dt, s, delayWindow, rateWindow); |
148 + | drfft0.FFT(); |
149 + | drfft0.searchPeak(); |
150 + | Float delay_out = drfft0.delay()(0, 1); |
151 + | Float rate_resn = 1.0f/(nPadfactor*nt*dt*1e9*f0); |
152 + | Float delay_resn = 1.0f/(nPadfactor*nchan*df); |
153 + | |
154 + | if (FRINGEJONES_TEST_VERBOSE) { |
155 + | cerr << boolalpha; |
156 + | cerr << "param = " << drfft0.param() << endl; |
157 + | cerr << "delay out = " << delay_out << endl; |
158 + | cerr << "delay in = " << delay << endl; |
159 + | cerr << "delay = " << delay << endl; |
160 + | cerr << "delay resolution = " << delay_resn << endl; |
161 + | cerr << "scaled error in delay = " << (delay_out - delay)/delay_resn << endl; |
162 + | cerr << "rate = " << drfft0.rate()(0, 1) << endl; |
163 + | cerr << "rate resolution = " << rate_resn << endl; |
164 + | cerr << "scaled error in rate = " << (drfft0.rate()(0, 1) - rate)/rate_resn << endl; |
165 + | } |
166 + | ASSERT_TRUE((delay_out < maxDelay) && (delay_out >= minDelay)); |
167 + | ASSERT_TRUE(allNearAbs(drfft0.rate()(0, 1), rate, 0.5*rate_resn)); |
168 + | } |
169 + | |
170 + | |
171 + | |
110 172 | class FringeJonesTest : public VisCalTestBase { |
111 173 | public: |
112 174 | // test values for solutions |
113 175 | Cube<Float> fpar; |
114 176 | |
115 177 | FringeJonesTest() : |
116 178 | VisCalTestBase(1, // nfield |
117 179 | 1, // nscan |
118 180 | 1, // nspw |
119 181 | 7, // nant |
123 185 | false), // unpolarized |
124 186 | // nPar, 1, {1 | nAntennas} |
125 187 | fpar(6, 1, VisCalTestBase::nAnt, 0.0f) // 6 pars per antenna |
126 188 | { |
127 189 | // Add FringeJonesTest specific init |
128 190 | // e.g., fill fpar with interesting values |
129 191 | for (Int i=1; i!=VisCalTestBase::nAnt; i++) { |
130 192 | // 1, 4 are delay. |
131 193 | fpar(1, 0, i) = 2.3; |
132 194 | fpar(4, 0, i) = -1.7; |
133 - | // 2 and 5 are rate. |
134 - | // |
135 - | // FIXME: We would like to use non-zero rates, but it seems |
136 - | // that the framework for corrupting the data might have some |
137 - | // missing pieces for FringeJones. |
138 - | // |
139 - | //fpar(2, 0, i) = -1.3e-9; |
140 - | // fpar(5, 0, i) = -1.7e-10; |
195 + | // Parameters 2 and 5 are rate. But VisCal::setMeta currently |
196 + | // only supports a single double for time, and this meta data is |
197 + | // used to generate all the Jones matrices for the FringeJones |
198 + | // (or any other VisCal subclass) calibrater; I don't feel |
199 + | // comfortable enough with the calibration stack to try to |
200 + | // intervene on this, so for now we can only test zero rates for |
201 + | // the FringeJones class proper. |
141 202 | fpar(2, 0, i) = 0.0; |
142 203 | fpar(5, 0, i) = 0.0; |
143 204 | } |
144 205 | // uncomment to see data shape summary from |
145 206 | //VisCalTestBase::summary("FringeJonesTest"); |
146 207 | } |
147 208 | }; |
148 209 | |
149 210 | |
150 - | |
151 211 | TEST_F(FringeJonesTest, FringeJonesApplyState) { |
152 212 | |
153 213 | FringeJones ff(VisCalTestBase::msmc); |
154 214 | ff.setApply(); |
155 215 | |
156 216 | ASSERT_EQ(VisCalEnum::JONES,ff.matrixType()); |
157 217 | ASSERT_EQ(VisCal::K,ff.type()); |
158 218 | ASSERT_EQ(String("Fringe Jones"),ff.typeName()); |
159 219 | //ASSERT_EQ(6,ff.nPar()); |
160 220 | ASSERT_FALSE(ff.freqDepPar()); |
175 235 | FringeJones FJapp(msmc); // "<noms>",nAnt,nSpw); |
176 236 | FJapp.setApply(); |
177 237 | // FJapp.setPrtlev(7); |
178 238 | |
179 239 | // Fill FJapp with actual parameters |
180 240 | for (Int ispw=0;ispw<nSpw;++ispw) { |
181 241 | FJapp.setMeta(0,1,0.0, |
182 242 | ispw,ssvp.freqs(ispw), |
183 243 | nChan); |
184 244 | FJapp.sizeApplyParCurrSpw(nChan); |
185 - | |
186 245 | // Enabled now phase model implemented... |
187 246 | FJapp.setApplyParCurrSpw(fpar,true,false); // don't invert |
188 247 | } |
189 248 | FringeJones FJsol(VisCalTestBase::msmc); |
190 249 | // FJsol.setPrtlev(7); |
191 250 | Record solvePar; |
192 251 | solvePar.define("table",String("test.Fringe")); // not used |
193 252 | solvePar.define("solint",String("inf")); |
194 253 | solvePar.define("combine",String("")); |
254 + | solvePar.define("globalsolve", true); |
255 + | solvePar.define("zerorates", true); |
195 256 | Vector<Int> refant(1,0); solvePar.define("refant",refant); |
257 + | Array<Double> delayWindow(IPosition(1, 2)); |
258 + | Array<Double> rateWindow(IPosition(1, 2)); |
259 + | delayWindow(IPosition(1, 0)) = -100.0; |
260 + | delayWindow(IPosition(1, 1)) = +100.0; |
261 + | rateWindow(IPosition(1, 0)) = -100.0; |
262 + | rateWindow(IPosition(1, 1)) = +100.0; |
263 + | solvePar.define("delaywindow", delayWindow); |
264 + | solvePar.define("ratewindow", rateWindow); |
265 + | |
196 266 | FJsol.setSolve(solvePar); |
197 267 | |
198 268 | SDBList sdbs; |
199 269 | Double reftime; |
200 270 | |
201 271 | Bool isFirst = true; |
202 272 | for (vi2.originChunks();vi2.moreChunks();vi2.nextChunk()) { |
203 273 | for (vi2.origin();vi2.more();vi2.next()) { |
204 274 | if (isFirst) { |
205 275 | reftime = vb2->time()(0); |
206 276 | isFirst = false; |
207 277 | } |
208 278 | Int ispw=vb2->spectralWindows()(0); |
209 279 | Int obsid(vb2->observationId()(0)); |
210 280 | Int scan(vb2->scan()(0)); |
211 - | Double timestamp(vb2->time()(0)); |
212 281 | Int fldid(vb2->fieldId()(0)); |
213 282 | Vector<Double> freqs(vb2->getFrequencies(0)); |
214 283 | Vector<Int> a1(vb2->antenna1()); |
215 284 | Vector<Int> a2(vb2->antenna2()); |
216 285 | |
217 286 | vb2->resetWeightsUsingSigma(); |
218 287 | vb2->setVisCubeCorrected(vb2->visCube()); |
219 288 | vb2->setFlagCube(vb2->flagCube()); |
220 289 | |
221 290 | // Corrupt with FJapp |
240 309 | FJsol.sizeSolveParCurrSpw(sdbs.nChannels()); |
241 310 | FJsol.selfSolveOne(sdbs); |
242 311 | |
243 312 | // Add comparison tests here |
244 313 | Matrix<Float> p = FJsol.solveRPar().nonDegenerate(); |
245 314 | Float delay1 = 2.3; |
246 315 | Float delay2 = -1.7; |
247 316 | Float rate1 = 0.0; |
248 317 | Float rate2 = 0.0; |
249 318 | |
319 + | if (FRINGEJONES_TEST_VERBOSE) { |
320 + | cerr << "delay1 results " << p(1,1) << endl; |
321 + | cerr << "delay2 results " << p(4,1) << endl; |
322 + | cerr << "rate1 results " << p(2,1) << endl; |
323 + | cerr << "rate2 results " << p(5,1) << endl; |
324 + | } |
250 325 | ASSERT_TRUE(allNearAbs(p(1, 1), delay1, 2e-2)); |
251 326 | ASSERT_TRUE(allNearAbs(p(4, 1), delay2, 2e-2)); |
252 327 | |
253 328 | ASSERT_TRUE(allNearAbs(p(2, 1), rate1, 1e-5)); |
254 329 | ASSERT_TRUE(allNearAbs(p(5, 1), rate2, 1e-5)); |
255 330 | |
256 331 | // cerr << "Parameters out: " << p << endl; |
257 332 | |
258 333 | } |