Line data Source code
1 : //# VisCalSolver.cc: Implementation of generic visibility solving
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 :
27 : #include <synthesis/MeasurementComponents/VisCalSolver.h>
28 :
29 : #include <msvis/MSVis/VisBuffer.h>
30 :
31 : #include <casacore/casa/Arrays/ArrayMath.h>
32 : #include <casacore/casa/Arrays/MaskArrMath.h>
33 : #include <casacore/casa/Arrays/ArrayLogical.h>
34 : #include <casacore/casa/Arrays/ArrayIter.h>
35 : //#include <scimath/Mathematics/MatrixMathLA.h>
36 : #include <casacore/casa/BasicSL/String.h>
37 : #include <casacore/casa/BasicMath/Math.h>
38 : #include <casacore/casa/Utilities/Assert.h>
39 : #include <casacore/casa/Exceptions/Error.h>
40 : #include <casacore/casa/OS/Memory.h>
41 : #include <casacore/casa/OS/Path.h>
42 :
43 : #include <sstream>
44 :
45 : #include <casacore/casa/Logging/LogMessage.h>
46 : #include <casacore/casa/Logging/LogSink.h>
47 :
48 : #define VCS_PRTLEV 0
49 :
50 : using namespace casacore;
51 : namespace casa { //# NAMESPACE CASA - BEGIN
52 :
53 :
54 : // **********************************************************
55 : // UVMod Implementations
56 : //
57 :
58 2 : VisCalSolver::VisCalSolver() :
59 : svb_(NULL),
60 : vbga_(NULL),
61 : ve_(NULL),
62 : svc_(NULL),
63 : nTotalPar_(0),
64 : R_(),dR_(),Rflg_(),
65 : maxIter_(50),
66 : chiSq_(0.0),
67 : chiSqV_(4,0.0),
68 : lastChiSq_(0.0),dChiSq_(0.0),
69 : sumWt_(0.0),nWt_(0),
70 : cvrgcount_(0),
71 : par_(), parOK_(), parErr_(), srcPar_(), lastCalPar_(), lastSrcPar_(),
72 : dpar_(), dcalpar_(), dsrcpar_(),
73 : grad_(),hess_(),
74 : lambda_(2.0),
75 : optstep_(true),
76 2 : prtlev_(VCS_PRTLEV)
77 : {
78 2 : if (prtlev()>0) cout << "VCS::VCS()" << endl;
79 2 : }
80 :
81 2 : VisCalSolver::~VisCalSolver()
82 : {
83 2 : if (prtlev()>0) cout << "VCS::~VCS()" << endl;
84 2 : }
85 :
86 0 : Bool VisCalSolver::solve(VisEquation& ve, SolvableVisCal& svc, VisBuffer& svb) {
87 :
88 0 : if (prtlev()>1) cout << "VCS::solve(,,)" << endl;
89 :
90 : /*
91 : LogSink logsink;
92 : {
93 : LogMessage message(LogOrigin("VisCalSolver", "solve"));
94 : ostringstream o; o<<"Beginning solve...";
95 : message.message(o);
96 : logsink.post(message);
97 : }
98 : */
99 : // Pointers to local ve,svc,svb
100 0 : ve_=&ve;
101 0 : svc_=&svc;
102 0 : svb_=&svb;
103 0 : vbga_=NULL;
104 :
105 : // Verify that VisEq has the correct svc:
106 : // TBD?
107 :
108 : // Initialize everything
109 0 : initSolve();
110 :
111 0 : Vector<Float> steplist(maxIter_+2,0.0);
112 0 : Vector<Float> rsteplist(maxIter_+2,0.0);
113 :
114 : // cout << "svb.modelVisCube() = " << phase(svb.modelVisCube())*180.0/C::pi << endl;
115 : // cout << "svb.modelVisCube() = " << amplitude(svb.modelVisCube()) << endl;
116 :
117 : // Verify VisBuffer validity for solving
118 : // (this sets parOK() on per-antenna basis (for focusChan)
119 : // based on data weights and baseline participation)
120 0 : Bool oktosolve = svc_->verifyForSolve(*svb_);
121 :
122 0 : if (oktosolve) {
123 :
124 0 : if (prtlev()>1) cout << "First guess:" << endl
125 0 : << "amp = " << amplitude(par()) << endl
126 0 : << "pha = " << phase(par())
127 0 : << endl;
128 :
129 : // Iterate solution
130 0 : Int iter(0);
131 0 : Bool done(false);
132 0 : while (!done) {
133 :
134 0 : if (prtlev()>2) cout << " Beginning iteration " << iter
135 0 : << "---------------------------------" << endl;
136 :
137 : // Differentiate the VB and get current Chi2
138 0 : differentiate();
139 0 : chiSquare();
140 0 : if (chiSq()==0.0) {
141 : // cout << "CHI2 IS SPURIOUSLY ZERO!*************************************" << endl;
142 : // cout << "R() = " << R() << endl;
143 : // cout << "wtmat = " << svb.weightMat() << endl;
144 : // cout << "flag = " << svb.flag() << endl;
145 : // cout << "sum(wtmat) = " << sum(wtmat) << endl;
146 0 : return false;
147 : }
148 :
149 0 : dChiSq() = chiSq()-lastChiSq();
150 :
151 : // cout << "chi2 = " << chiSq() << " " << dChiSq() << " " << dChiSq()/chiSq() << endl;
152 :
153 : // Continuue if we haven't converged
154 0 : if (!converged()) {
155 :
156 0 : if (dChiSq()<=0.0) {
157 : // last step was good...
158 0 : lastChiSq()=chiSq();
159 :
160 : // so accumulate new grad/hess...
161 0 : accGradHess();
162 :
163 : //...and adjust lambda downward
164 : // lambda()/=2.0;
165 : // lambda()=0.8;
166 0 : lambda()=1.0;
167 : }
168 : else {
169 : // cout << "reverting..." << chiSq() << " " << dChiSq() << " (" << iter << ")" << endl;
170 : // last step was bad, revert to previous
171 0 : revert();
172 : //...with a larger lambda
173 : // lambda()*=4.0;
174 0 : lambda()=1.0;
175 : }
176 :
177 : // Solve for the parameter step
178 0 : solveGradHess();
179 :
180 : // Remember curr pars
181 0 : lastCalPar()=par();
182 0 : if (svc_->solvePol())
183 0 : lastSrcPar()=srcPar();
184 :
185 : // Refine the step size by exploring chi2 in the
186 : // gradient direction
187 0 : if (optstep_) // && cvrgcount_>=3)
188 0 : optStepSize();
189 :
190 : // Update current parameters (saves a copy of them)
191 0 : updatePar();
192 :
193 : // cout << "srcPar() = " << srcPar() << endl;
194 :
195 0 : steplist(iter)=max(amplitude(dCalPar()));
196 0 : rsteplist(iter)=max(amplitude(dCalPar())/amplitude(par()));
197 :
198 : }
199 : else {
200 : // Convergence means we're done!
201 0 : done=true;
202 :
203 0 : if (prtlev()>0) {
204 0 : cout << "Iterations =" << iter << endl;
205 0 : cout << "par()=" << par() << endl;
206 0 : cout << "srcPar()=" << srcPar() << endl;
207 : }
208 :
209 :
210 : /*
211 : cout << " good pars=" << ntrue(parOK())
212 : << " iterations=" << iter << endl
213 : << " steps=" << steplist(IPosition(1,0),IPosition(1,iter))
214 : << endl
215 : << " rsteps=" << rsteplist(IPosition(1,0),IPosition(1,iter))
216 : << endl;
217 : */
218 :
219 : // Get parameter errors:
220 0 : accGradHess();
221 0 : getErrors();
222 :
223 : // Return, signaling success if at least 1 good solution
224 0 : return (ntrue(parOK())>0);
225 :
226 : }
227 :
228 : // Escape iteration loop via iteration limit
229 0 : if (iter==maxIter()) {
230 0 : cout << "Reached iteration limit: " << iter << " iterations. " << endl;
231 : // cout << " good pars = " << ntrue(parOK())
232 : // << " steps = " << steplist
233 : // << endl;
234 0 : done=true;
235 : }
236 :
237 : // Advance iteration counter
238 0 : iter++;
239 : }
240 :
241 : }
242 : else {
243 0 : cout << " Insufficient unflagged antennas to proceed with this solve." << endl;
244 : }
245 :
246 0 : return false;
247 :
248 : }
249 :
250 : // New VisBuffGroupAcc version
251 0 : Bool VisCalSolver::solve(VisEquation& ve, SolvableVisCal& svc, VisBuffGroupAcc& vbga) {
252 :
253 0 : if (prtlev()>1) cout << "VCS::solve(,,VBGA)" << endl;
254 :
255 : /*
256 : LogSink logsink;
257 : {
258 : LogMessage message(LogOrigin("VisCalSolver", "solve"));
259 : ostringstream o; o<<"Beginning solve...";
260 : message.message(o);
261 : logsink.post(message);
262 : }
263 : */
264 : // Pointers to local ve,svc,svb
265 0 : ve_=&ve;
266 0 : svc_=&svc;
267 0 : svb_=NULL;
268 0 : vbga_=&vbga;
269 :
270 : // Verify that VisEq has the correct svc:
271 : // TBD?
272 :
273 : // Initialize everything
274 0 : initSolve();
275 :
276 0 : Vector<Float> steplist(maxIter_+2,0.0);
277 0 : Vector<Float> rsteplist(maxIter_+2,0.0);
278 :
279 : // Verify Data's validity for solve w.r.t. baselines available
280 : // (this sets parOK() on per-antenna basis (for focusChan)
281 : // based on data weights and baseline participation)
282 0 : Bool oktosolve = svc_->verifyConstraints(*vbga_);
283 :
284 0 : if (oktosolve) {
285 :
286 0 : if (prtlev()>1) cout << "First guess:" << endl
287 0 : << "amp = " << amplitude(par()) << endl
288 0 : << "pha = " << phase(par())
289 0 : << endl;
290 :
291 : // Iterate solution
292 0 : Int iter(0);
293 0 : Bool done(false);
294 0 : while (!done) {
295 :
296 0 : if (prtlev()>2) cout << " Beginning iteration " << iter
297 0 : << "---------------------------------" << endl;
298 :
299 : // Differentiate the VB and get current Chi2
300 0 : differentiate2();
301 0 : chiSquare2();
302 0 : if (chiSq()==0.0) {
303 : // cout << "CHI2 IS SPURIOUSLY ZERO!*************************************" << endl;
304 : // cout << "R() = " << R() << endl;
305 : // cout << "wtmat = " << svb.weightMat() << endl;
306 : // cout << "flag = " << svb.flag() << endl;
307 : // cout << "sum(wtmat) = " << sum(wtmat) << endl;
308 0 : return false;
309 : }
310 :
311 0 : dChiSq() = chiSq()-lastChiSq();
312 :
313 : // cout << "chi2 = " << chiSq() << " " << dChiSq() << " " << dChiSq()/chiSq() << endl;
314 :
315 : // Continuue if we haven't converged
316 0 : if (!converged()) {
317 :
318 0 : if (dChiSq()<=0.0) {
319 : // last step was good...
320 0 : lastChiSq()=chiSq();
321 :
322 : // so accumulate new grad/hess...
323 0 : accGradHess2();
324 :
325 : //...and adjust lambda downward
326 : // lambda()/=2.0;
327 : // lambda()=0.8;
328 0 : lambda()=1.0;
329 : }
330 : else {
331 : // cout << "reverting..." << chiSq() << " " << dChiSq() << " (" << iter << ")" << endl;
332 : // last step was bad, revert to previous
333 0 : revert();
334 : //...with a larger lambda
335 : // lambda()*=4.0;
336 0 : lambda()=1.0;
337 : }
338 :
339 : // Solve for the parameter step
340 0 : solveGradHess();
341 :
342 : // Remember curr pars
343 0 : lastCalPar()=par();
344 0 : if (svc_->solvePol())
345 0 : lastSrcPar()=srcPar();
346 :
347 : // Refine the step size by exploring chi2 in the
348 : // gradient direction
349 0 : if (optstep_) // && cvrgcount_>=3)
350 0 : optStepSize2();
351 :
352 : // Update current parameters (saves a copy of them)
353 0 : updatePar();
354 :
355 : // cout << "srcPar() = " << srcPar() << endl;
356 :
357 0 : steplist(iter)=max(amplitude(dCalPar()));
358 0 : rsteplist(iter)=max(amplitude(dCalPar())/amplitude(par()));
359 :
360 : }
361 : else {
362 : // Convergence means we're done!
363 0 : done=true;
364 :
365 0 : if (prtlev()>0) {
366 0 : cout << "Iterations =" << iter << endl;
367 0 : cout << "par()=" << par() << endl;
368 0 : cout << "srcPar()=" << srcPar() << endl;
369 : }
370 :
371 :
372 : /*
373 : cout << " good pars=" << ntrue(parOK())
374 : << " iterations=" << iter << endl
375 : << " steps=" << steplist(IPosition(1,0),IPosition(1,iter))
376 : << endl
377 : << " rsteps=" << rsteplist(IPosition(1,0),IPosition(1,iter))
378 : << endl;
379 : */
380 :
381 : // Get parameter errors:
382 0 : accGradHess2();
383 0 : getErrors();
384 :
385 : // Return, signaling success if at least 1 good solution
386 0 : return (ntrue(parOK())>0);
387 :
388 : }
389 :
390 : // Escape iteration loop via iteration limit
391 0 : if (iter==maxIter()) {
392 0 : cout << "Reached iteration limit: " << iter << " iterations. " << endl;
393 : // cout << " good pars = " << ntrue(parOK())
394 : // << " steps = " << steplist
395 : // << endl;
396 0 : done=true;
397 : }
398 :
399 : // Advance iteration counter
400 0 : iter++;
401 : }
402 :
403 : }
404 : else {
405 0 : cout << " Insufficient unflagged antennas to proceed with this solve." << endl;
406 : }
407 :
408 0 : return false;
409 :
410 : }
411 :
412 0 : void VisCalSolver::initSolve() {
413 :
414 0 : if (prtlev()>2) cout << " VCS::initSolve()" << endl;
415 :
416 : // Get total number of cal parameters from svc info
417 0 : nCalPar()=svc().nTotalPar();
418 :
419 : // solvePol() tells us how many source pol params we are solving for
420 0 : nSrcPar()=svc().solvePol();
421 :
422 : // the total number of parameters
423 0 : nTotalPar()=nCalPar()+nSrcPar();
424 :
425 0 : if (prtlev()>2)
426 0 : cout << " Total parameters in solve: " << nTotalPar() << endl;
427 :
428 : // Chi2 and weights
429 0 : chiSq()=0.0;
430 0 : lastChiSq()=DBL_MAX;
431 0 : dChiSq()=0.0;
432 :
433 0 : sumWt()=0.0;
434 0 : nWt()=0;
435 :
436 : // Link up svc's internal pars with local reference
437 : // (only if shape is correct)
438 :
439 0 : if (svc().solveCPar().nelements()==uInt(nCalPar())) {
440 0 : par().reference(svc().solveCPar().reform(IPosition(1,nCalPar())));
441 0 : parOK().reference(svc().solveParOK().reform(IPosition(1,nCalPar())));
442 0 : parErr().reference(svc().solveParErr().reform(IPosition(1,nCalPar())));
443 0 : if (svc().solvePol())
444 0 : srcPar().reference(svc().srcPolPar());
445 : }
446 : else
447 0 : throw(AipsError("Solver and SVC cannot synchronize parameters."));
448 :
449 : // Pars
450 :
451 0 : dpar().resize(nTotalPar());
452 0 : dpar()=0.0;
453 :
454 0 : lastCalPar().resize(nCalPar());
455 0 : dCalPar().reference(dpar()(IPosition(1,0),IPosition(1,nCalPar()-1)));
456 :
457 0 : if (svc().solvePol()) {
458 0 : lastSrcPar().resize(nSrcPar());
459 0 : dSrcPar().reference(dpar()(IPosition(1,nCalPar()),IPosition(1,nTotalPar()-1)));
460 : }
461 : else {
462 0 : lastSrcPar().resize();
463 0 : dSrcPar().resize();
464 : }
465 :
466 : // Gradient and Hessian
467 0 : grad().resize(nTotalPar());
468 0 : grad()=0.0;
469 :
470 0 : hess().resize(nTotalPar());
471 0 : hess()=0.0;
472 :
473 : // Levenberg-Marquardt factor
474 0 : lambda()=2.0;
475 :
476 : // Convergence anticipation
477 0 : cvrgcount_=0;
478 :
479 0 : }
480 :
481 0 : void VisCalSolver::residualate() {
482 :
483 0 : if (prtlev()>2) cout << " VCS::residualate()" << endl;
484 :
485 : // Delegate to VisEquation
486 : // ve().residuals(svb(),R(),Rflag());
487 :
488 : // For now, just use ve.diffResid, until we have
489 : // implemented focuschan-aware trial corrupt in SVC
490 : // (this will hurt performance a bit)
491 0 : ve().diffResiduals(svb(),R(),dR(),Rflg());
492 :
493 0 : }
494 :
495 0 : void VisCalSolver::residualate2() {
496 :
497 0 : if (prtlev()>2) cout << " VCS::residualate()" << endl;
498 :
499 : // Delegate to VisEquation
500 : // ve().residuals(svb(),R(),Rflag());
501 :
502 : // For now, just use ve.diffResid, until we have
503 : // implemented focuschan-aware trial corrupt in SVC
504 : // (this will hurt performance a bit)
505 :
506 0 : for (Int ibuf=0;ibuf<vbga().nBuf();++ibuf)
507 0 : ve().diffResiduals(vbga()(ibuf));
508 0 : }
509 :
510 0 : void VisCalSolver::differentiate() {
511 :
512 0 : if (prtlev()>2) cout << " VCS::differentiate()" << endl;
513 :
514 : // Delegate to VisEquation
515 0 : ve().diffResiduals(svb(),R(),dR(),Rflg());
516 :
517 :
518 0 : if (svc().solvePol()) {
519 :
520 : // Differentiate w.r.t source
521 0 : svc().diffSrc(svb(),dSrc());
522 : }
523 :
524 0 : if (prtlev()>6) { // R, dR
525 0 : cout << " R= " << R() << endl;
526 0 : cout << " dR=" << dR() << endl;
527 : }
528 :
529 0 : }
530 :
531 0 : void VisCalSolver::differentiate2() {
532 :
533 0 : if (prtlev()>2) cout << " VCS::differentiate(VBGA version)" << endl;
534 :
535 : // Delegate to VisEquation
536 0 : for (Int ibuf=0;ibuf<vbga().nBuf();++ibuf)
537 0 : ve().diffResiduals(vbga()(ibuf));
538 :
539 0 : if (svc().solvePol()) {
540 : // throw(AipsError("solvePol not yet sorted w.r.t. VBGA."));
541 : // Differentiate w.r.t source
542 0 : svc().diffSrc(vbga()(0),dSrc());
543 : }
544 :
545 0 : }
546 :
547 0 : void VisCalSolver::chiSquare() {
548 :
549 0 : if (prtlev()>2) cout << " VCS::chiSquare()" << endl;
550 :
551 : // NB: Assumes R() is up-to-date
552 :
553 : // TBD: Review correctness of summing weights
554 : // inside the channel loop?
555 :
556 : // TBD: per-ant/bln chiSq?
557 :
558 : // cout << "VCS::chiSquare: svc().focusChan() = " << svc().focusChan() << endl;
559 :
560 0 : chiSq()=0.0;
561 0 : chiSqV()=0.0;
562 0 : sumWt()=0.0;
563 0 : nWt()=0;
564 :
565 : // Shapes for iteration
566 0 : IPosition shR(R().shape());
567 0 : Int nCorr=shR(0);
568 0 : Int nChan=shR(1);
569 0 : Int nRow=shR(2);
570 :
571 0 : ArrayIterator<Bool> flag(svb().flag(),1); // fl(chan) by (row)
572 0 : ArrayIterator<Float> wtMat(svb().weightMat(),1); // wt(corr) by (row)
573 0 : ArrayIterator<Complex> Rit(R(),1); // R(corr) by (chan,row)
574 :
575 0 : Bool* flR = svb().flagRow().data();
576 : Bool* fl;
577 : Float* wt;
578 : Complex *Rp;
579 :
580 0 : for (Int irow=0;irow<nRow;++irow) {
581 0 : if (!*flR) {
582 : // This row's wt(corr), flag(chan)
583 0 : wt = wtMat.array().data();
584 : // wt[1]=wt[2]=0.0;
585 0 : fl = flag.array().data() + svc().focusChan();
586 : // Register R for this row, 0th channel
587 0 : for (Int ich=0;ich<nChan;++ich) {
588 0 : if (!*fl) {
589 :
590 0 : Rp = Rit.array().data();
591 0 : for (Int icorr=0;icorr<nCorr;++icorr) {
592 :
593 : /*
594 : if (svb().weightMat()(icorr,irow)>0.0) {
595 : cout << irow << " " << icorr << " "
596 : << svb().weightMat()(icorr,irow) << " "
597 : << *wt << " "
598 : << R()(icorr,0,irow) << " "
599 : << *Rp << " "
600 : << &R()(icorr,0,irow) << " "
601 : << Rp << " "
602 : << Rp-&R()(icorr,0,irow) << " "
603 : << &R()(0,0,irow)-&R()(3,0,irow) << " "
604 : << Rit.array().shape()
605 : << endl;
606 : }
607 : */
608 0 : chiSq()+=Double( (*wt)*real((*Rp)*conj(*Rp)) );
609 : // chiSq()+=Double( real((*Rp)*conj(*Rp)) );
610 0 : chiSqV()(icorr)+=Double( (*wt)*real((*Rp)*conj(*Rp)) );
611 0 : sumWt()+=Double(*wt); // for each channel?!
612 :
613 0 : if (*wt>0.0) nWt()++;
614 :
615 : //Advance to next corr
616 0 : ++wt;
617 0 : ++Rp;
618 : }
619 : // Use same wt(corr) vectors for each channel
620 0 : wt-=nCorr;
621 : }
622 : // Advance to next channel
623 0 : Rit.next();
624 0 : ++fl;
625 : }
626 : }
627 : else
628 : // Advance over flagged row!
629 0 : for (Int ich=0;ich<nChan;++ich) Rit.next();
630 :
631 : // Advance to next row
632 0 : ++flR;
633 0 : flag.next();
634 0 : wtMat.next();
635 : }
636 :
637 :
638 0 : }
639 :
640 0 : void VisCalSolver::chiSquare2() {
641 :
642 0 : if (prtlev()>2) cout << " VCS::chiSquare(CVB version)" << endl;
643 :
644 : // NB: Assumes R() is up-to-date
645 :
646 : // TBD: Review correctness of summing weights
647 : // inside the channel loop?
648 :
649 : // TBD: per-ant/bln chiSq?
650 :
651 : // cout << "VCS::chiSquare: svc().focusChan() = " << svc().focusChan() << endl;
652 :
653 0 : chiSq()=0.0;
654 0 : chiSqV()=0.0;
655 0 : sumWt()=0.0;
656 0 : nWt()=0;
657 :
658 : // Loop over CVBs
659 0 : for (Int ibuf=0;ibuf<vbga().nBuf();++ibuf) {
660 :
661 : // Focus on one CVB at a time
662 0 : CalVisBuffer& cvb(vbga()(ibuf));
663 :
664 0 : R().reference(cvb.residuals());
665 :
666 : // Shapes for iteration
667 0 : IPosition shR(R().shape());
668 0 : Int nCorr=shR(0);
669 0 : Int nChan=shR(1);
670 0 : Int nRow=shR(2);
671 :
672 0 : ArrayIterator<Bool> flag(cvb.residFlag(),1); // fl(chan) by (row)
673 0 : ArrayIterator<Float> wtMat(cvb.weightMat(),1); // wt(corr) by (row)
674 0 : ArrayIterator<Complex> Rit(R(),1); // R(corr) by (chan,row)
675 :
676 0 : Bool* flR = cvb.flagRow().data();
677 : Bool* fl;
678 : Float* wt;
679 : Complex *Rp;
680 :
681 0 : for (Int irow=0;irow<nRow;++irow) {
682 0 : if (!*flR) {
683 : // This row's wt(corr), flag(chan)
684 0 : wt = wtMat.array().data();
685 0 : fl = flag.array().data();
686 : // Register R for this row, 0th channel
687 0 : for (Int ich=0;ich<nChan;++ich) {
688 0 : if (!*fl) {
689 :
690 0 : Rp = Rit.array().data();
691 0 : for (Int icorr=0;icorr<nCorr;++icorr) {
692 :
693 0 : chiSq()+=Double( (*wt)*real((*Rp)*conj(*Rp)) );
694 0 : chiSqV()(icorr)+=Double( (*wt)*real((*Rp)*conj(*Rp)) );
695 0 : sumWt()+=Double(*wt); // for each channel?!
696 :
697 0 : if (*wt>0.0) nWt()++;
698 :
699 : //Advance to next corr
700 0 : ++wt;
701 0 : ++Rp;
702 : }
703 : // Use same wt(corr) vectors for each channel
704 0 : wt-=nCorr;
705 : }
706 : // Advance to next channel
707 0 : Rit.next();
708 0 : ++fl;
709 : }
710 : }
711 : else
712 : // Advance over flagged row!
713 0 : for (Int ich=0;ich<nChan;++ich) Rit.next();
714 :
715 : // Advance to next row
716 0 : ++flR;
717 0 : flag.next();
718 0 : wtMat.next();
719 : }
720 : } // ibuf
721 0 : }
722 :
723 :
724 0 : Bool VisCalSolver::converged() {
725 :
726 0 : if (prtlev()>2) cout << " VCS::converged()" << endl;
727 :
728 : // Change in chi2
729 0 : dChiSq() = chiSq()-lastChiSq();
730 0 : Float fChiSq(dChiSq()/chiSq());
731 :
732 : // if (prtlev()>0)
733 : // cout << "chiSq: " << chiSq() << " " << chiSqV() << " " << lambda() << endl;
734 :
735 : // Consider convergence if chi2 decreases...
736 : // if (dChiSq()<=0.0) {
737 0 : if (fChiSq<=0.001) {
738 :
739 : // ...and the change is small:
740 0 : if (abs(dChiSq()) < 0.1*chiSq()) {
741 0 : ++cvrgcount_;
742 :
743 : // if (cvrgcount_==2) lambda()=2.0;
744 :
745 : // Four such steps we believe we have converged!
746 : // if (cvrgcount_>3)
747 : //if (cvrgcount_>10)
748 0 : if (cvrgcount_>5)
749 0 : return true;
750 : }
751 :
752 0 : if (prtlev()>0)
753 0 : cout << " Good: chiSq=" << chiSq()
754 0 : << " dChiSq=" << dChiSq()
755 0 : << " fChiSq=" << dChiSq()/chiSq()
756 0 : << " cvrgcnt=" << cvrgcount_
757 0 : << " lambda=" << lambda()
758 0 : << endl;
759 :
760 : }
761 : else {
762 : // (chi2 increased)
763 :
764 : // If a large increase, don't anticipate yet
765 0 : if (abs(dChiSq()) > 0.1*chiSq())
766 0 : cvrgcount_=0;
767 : else {
768 : // anticipate a little less if upward change is small
769 : // TBD: is this right?
770 0 : --cvrgcount_;
771 0 : cvrgcount_=max(cvrgcount_,0); // never less than zero
772 : }
773 :
774 0 : if (prtlev()>0)
775 0 : cout << " Bad: chiSq=" << chiSq()
776 0 : << " dChiSq=" << dChiSq()
777 0 : << " fChiSq=" << dChiSq()/chiSq()
778 0 : << " cvrgcnt=" << cvrgcount_
779 0 : << " lambda=" << lambda()
780 0 : << endl;
781 :
782 :
783 : }
784 :
785 : // Not yet converged
786 0 : return false;
787 :
788 : }
789 :
790 :
791 0 : void VisCalSolver::accGradHess() {
792 :
793 0 : if (prtlev()>2) cout << " VCS::accGradHess()" << endl;
794 :
795 0 : grad()=0.0;
796 0 : hess()=0.0;
797 :
798 0 : IPosition dRip(dR().shape());
799 :
800 0 : Int nRow(dRip(3));
801 0 : Int nChan(dRip(2));
802 0 : Int nPar(dRip(1)); // pars per antenna
803 0 : Int nCorr(dRip(0));
804 :
805 0 : ArrayIterator<Bool> flag(svb().flag(),1); // fl(chan) by (row)
806 0 : ArrayIterator<Float> wtMat(svb().weightMat(),1); // wt(corr) by (row)
807 :
808 0 : ArrayIterator<Complex> Rit(R(),1); // R(corr) by (chan,row)
809 :
810 0 : Array<Complex> dR0, dR1;
811 : {
812 0 : ArrayIterator<Complex> dRit(dR(),4); // dR(corr,par,chan,row) by (ant)
813 0 : dR0.reference(dRit.array());
814 0 : dRit.next();
815 0 : dR1.reference(dRit.array());
816 : }
817 0 : ArrayIterator<Complex> dR0it(dR0,2); // dR0(corr,par) by (chan,row)
818 0 : ArrayIterator<Complex> dR1it(dR1,2); // dR1(corr,par) by (chan,row)
819 :
820 0 : Bool* flR = svb().flagRow().data();
821 : Bool* fl;
822 : Float* wt;
823 0 : Int* a1 = svb().antenna1().data();
824 0 : Int* a2 = svb().antenna2().data();
825 : Complex *Rp;
826 : Complex *dR0p;
827 : Complex *dR1p;
828 : DComplex *G1;
829 : DComplex *G2;
830 : Double *H1;
831 : Double *H2;
832 :
833 : // cout << "flag.array().data() = " << flag.array().data() << endl;
834 : // cout << "VCS::accGradHess: svc().focusChan() = " << svc().focusChan() << endl;
835 :
836 0 : for (Int irow=0;irow<nRow;++irow) {
837 0 : if (!*flR) {
838 : // Register grad, hess for ants in this baseline
839 0 : Int p1((*a1)*nPar);
840 0 : Int p2((*a2)*nPar);
841 0 : G1 = grad().data() + p1;
842 0 : G2 = grad().data() + p2;
843 0 : H1 = hess().data() + p1;
844 0 : H2 = hess().data() + p2;
845 : // This row's wt(corr), flag(chan)
846 0 : wt = wtMat.array().data();
847 0 : fl = flag.array().data() + svc().focusChan();
848 0 : for (Int ich=0;ich<nChan;++ich) {
849 0 : if (!*fl) {
850 :
851 : // Do source bits, if necessary
852 0 : if (svc().solvePol()) {
853 0 : for (Int icorr=0;icorr<nCorr;++icorr) {
854 0 : Float swt=Vector<Float>(wtMat.array())(icorr);
855 :
856 0 : for (Int ispar=0;ispar<nSrcPar();++ispar) {
857 0 : grad()(nCalPar()+ispar) += DComplex((swt)*2.0*real(dSrc()(IPosition(4,icorr,ich,irow,ispar))*
858 0 : conj(R()(icorr,ich,irow))));
859 0 : hess()(nCalPar()+ispar) += Double((swt)*2.0*real(dSrc()(IPosition(4,icorr,ich,irow,ispar))*
860 0 : conj(dSrc()(IPosition(4,icorr,ich,irow,ispar)))));
861 : }
862 : }
863 : }
864 :
865 : // Register R,dR for this channel, row
866 0 : Rp = Rit.array().data();
867 0 : dR0p = dR0it.array().data();
868 0 : dR1p = dR1it.array().data();
869 :
870 : /*
871 : cout << "irow=" << irow << endl;
872 : cout << "R=" << Rit.array() << endl;
873 : cout << "dR0=" << dR0it.array()
874 : << "dR1=" << dR1it.array()
875 : << endl;
876 : */
877 :
878 0 : for (int ip=0;ip<nPar;++ip) {
879 :
880 0 : for (Int icorr=0;icorr<nCorr;++icorr) {
881 :
882 0 : (*G1) += DComplex( (*wt)*((*Rp) *conj(*dR0p)) );
883 0 : (*G2) += DComplex( (*wt)*((*dR1p)*conj(*Rp)) );
884 0 : (*H1) += Double( (*wt)*real((*dR0p)*conj(*dR0p)) );
885 0 : (*H2) += Double( (*wt)*real((*dR1p)*conj(*dR1p)) );
886 : /*
887 : (*G1) += DComplex( (*Rp) *conj(*dR0p) );
888 : (*G2) += DComplex( (*dR1p)*conj(*Rp) );
889 : (*H1) += Double( real((*dR0p)*conj(*dR0p)) );
890 : (*H2) += Double( real((*dR1p)*conj(*dR1p)) );
891 : */
892 :
893 :
894 : //Advance to next corr
895 0 : ++wt;
896 0 : ++Rp;
897 0 : ++dR0p;++dR1p;
898 : }
899 : // Advance to next par
900 0 : ++G1; ++G2;
901 0 : ++H1; ++H2;
902 : // Use same Rp(corr), wt(corr) vectors for each parameter
903 0 : Rp-=nCorr;
904 0 : wt-=nCorr;
905 : }
906 : // Accumulate to same grad(par) & hess(par) for each channel
907 0 : G1-=nPar; G2-=nPar;
908 0 : H1-=nPar; H2-=nPar;
909 : }
910 : // Advance to next channel
911 0 : ++fl;
912 0 : Rit.next();
913 0 : dR0it.next();
914 0 : dR1it.next();
915 : }
916 : } // !*flgR
917 : else {
918 : // Advance over flagged row
919 0 : for (Int ich=0;ich<nChan;++ich) {
920 0 : Rit.next();
921 0 : dR0it.next();
922 0 : dR1it.next();
923 : }
924 : }
925 : // Advance to next row
926 0 : ++flR;
927 0 : ++a1;++a2;
928 0 : flag.next();
929 0 : wtMat.next();
930 : }
931 :
932 0 : if (prtlev()>4) { // grad, hess
933 0 : cout << " grad= " << grad() << endl;
934 0 : cout << " hess= " << hess() << endl;
935 : }
936 :
937 :
938 : /*
939 : // This is the slower way with array indexing
940 : for (Int irow=0;irow<nRow;++irow,++flR,++a1,++a2) {
941 : if (!*flR) {
942 : par1=(*a1)*nPar;
943 : par2=(*a2)*nPar;
944 : for (Int ich=0;ich<nChan;++ich,++fl) {
945 : if (!*fl) {
946 : for (int ip=0;ip<nPar;++ip) {
947 : for (Int icorr=0;icorr<nCorr;++icorr) {
948 : grad(par1+ip) += ( Double(wt(icorr,irow))*
949 : Double( R()(icorr,ich,irow )*
950 : conj(dR()(icorr,ip,ich,irow,0)) ) );
951 : grad(par2+ip) += ( Double(wt(icorr,irow))*
952 : Double( dR()(icorr,ip,ich,irow,1)*
953 : conj(R()(icorr,ich,irow )) ) );
954 :
955 : hess(par1+ip) += ( Double(wt(icorr,irow))*
956 : Double( dR()(icorr,ip,ich,irow,0)*
957 : conj(dR()(icorr,ip,ich,irow,0)) ) );
958 : hess(par2+ip) += ( Double(wt(icorr,irow))*
959 : Double( dR()(icorr,ip,ich,irow,1)*
960 : conj(dR()(icorr,ip,ich,irow,1)) ) );
961 : }
962 : }
963 : }
964 : }
965 : }
966 : else {
967 : // Advance over flagged row
968 : fl+=nChan;
969 : }
970 :
971 : }
972 : */
973 :
974 :
975 0 : }
976 :
977 :
978 0 : void VisCalSolver::accGradHess2() {
979 :
980 0 : if (prtlev()>2) cout << " VCS::accGradHess(CVB version)" << endl;
981 :
982 0 : grad()=0.0;
983 0 : hess()=0.0;
984 :
985 0 : for (Int ibuf=0;ibuf<vbga().nBuf();++ibuf) {
986 :
987 0 : CalVisBuffer& cvb(vbga()(ibuf));
988 :
989 0 : R().reference(cvb.residuals());
990 0 : dR().reference(cvb.diffResiduals());
991 :
992 :
993 0 : IPosition dRip(dR().shape());
994 :
995 0 : Int nRow(dRip(3));
996 0 : Int nChan(dRip(2));
997 0 : Int nPar(dRip(1)); // pars per antenna
998 0 : Int nCorr(dRip(0));
999 :
1000 0 : ArrayIterator<Bool> flag(cvb.residFlag(),1); // fl(chan) by (row)
1001 0 : ArrayIterator<Float> wtMat(cvb.weightMat(),1); // wt(corr) by (row)
1002 :
1003 0 : ArrayIterator<Complex> Rit(R(),1); // R(corr) by (chan,row)
1004 :
1005 0 : Array<Complex> dR0, dR1;
1006 : {
1007 0 : ArrayIterator<Complex> dRit(dR(),4); // dR(corr,par,chan,row) by (ant)
1008 0 : dR0.reference(dRit.array());
1009 0 : dRit.next();
1010 0 : dR1.reference(dRit.array());
1011 : }
1012 0 : ArrayIterator<Complex> dR0it(dR0,2); // dR0(corr,par) by (chan,row)
1013 0 : ArrayIterator<Complex> dR1it(dR1,2); // dR1(corr,par) by (chan,row)
1014 :
1015 0 : Bool* flR = cvb.flagRow().data();
1016 : Bool* fl;
1017 : Float* wt;
1018 0 : Int* a1 = cvb.antenna1().data();
1019 0 : Int* a2 = cvb.antenna2().data();
1020 : Complex *Rp;
1021 : Complex *dR0p;
1022 : Complex *dR1p;
1023 : DComplex *G1;
1024 : DComplex *G2;
1025 : Double *H1;
1026 : Double *H2;
1027 :
1028 0 : for (Int irow=0;irow<nRow;++irow) {
1029 0 : if (!*flR) {
1030 : // Register grad, hess for ants in this baseline
1031 0 : Int p1((*a1)*nPar);
1032 0 : Int p2((*a2)*nPar);
1033 0 : G1 = grad().data() + p1;
1034 0 : G2 = grad().data() + p2;
1035 0 : H1 = hess().data() + p1;
1036 0 : H2 = hess().data() + p2;
1037 : // This row's wt(corr), flag(chan)
1038 0 : wt = wtMat.array().data();
1039 0 : fl = flag.array().data();
1040 0 : for (Int ich=0;ich<nChan;++ich) {
1041 0 : if (!*fl) {
1042 :
1043 : // Do source bits, if necessary
1044 0 : if (svc().solvePol()) {
1045 0 : for (Int icorr=0;icorr<nCorr;++icorr) {
1046 0 : Float swt=Vector<Float>(wtMat.array())(icorr);
1047 :
1048 0 : for (Int ispar=0;ispar<nSrcPar();++ispar) {
1049 0 : grad()(nCalPar()+ispar) += DComplex((swt)*2.0*real(dSrc()(IPosition(4,icorr,ich,irow,ispar))*
1050 0 : conj(R()(icorr,ich,irow))));
1051 0 : hess()(nCalPar()+ispar) += Double((swt)*2.0*real(dSrc()(IPosition(4,icorr,ich,irow,ispar))*
1052 0 : conj(dSrc()(IPosition(4,icorr,ich,irow,ispar)))));
1053 : }
1054 : }
1055 : }
1056 :
1057 : // Register R,dR for this channel, row
1058 0 : Rp = Rit.array().data();
1059 0 : dR0p = dR0it.array().data();
1060 0 : dR1p = dR1it.array().data();
1061 :
1062 0 : for (int ip=0;ip<nPar;++ip) {
1063 :
1064 0 : for (Int icorr=0;icorr<nCorr;++icorr) {
1065 :
1066 0 : (*G1) += DComplex( (*wt)*((*Rp) *conj(*dR0p)) );
1067 0 : (*G2) += DComplex( (*wt)*((*dR1p)*conj(*Rp)) );
1068 0 : (*H1) += Double( (*wt)*real((*dR0p)*conj(*dR0p)) );
1069 0 : (*H2) += Double( (*wt)*real((*dR1p)*conj(*dR1p)) );
1070 :
1071 : //Advance to next corr
1072 0 : ++wt;
1073 0 : ++Rp;
1074 0 : ++dR0p;++dR1p;
1075 : }
1076 : // Advance to next par
1077 0 : ++G1; ++G2;
1078 0 : ++H1; ++H2;
1079 : // Use same Rp(corr), wt(corr) vectors for each parameter
1080 0 : Rp-=nCorr;
1081 0 : wt-=nCorr;
1082 : }
1083 : // Accumulate to same grad(par) & hess(par) for each channel
1084 0 : G1-=nPar; G2-=nPar;
1085 0 : H1-=nPar; H2-=nPar;
1086 : }
1087 : // Advance to next channel
1088 0 : ++fl;
1089 0 : Rit.next();
1090 0 : dR0it.next();
1091 0 : dR1it.next();
1092 : }
1093 : } // !*flgR
1094 : else {
1095 : // Advance over flagged row
1096 0 : for (Int ich=0;ich<nChan;++ich) {
1097 0 : Rit.next();
1098 0 : dR0it.next();
1099 0 : dR1it.next();
1100 : }
1101 : }
1102 : // Advance to next row
1103 0 : ++flR;
1104 0 : ++a1;++a2;
1105 0 : flag.next();
1106 0 : wtMat.next();
1107 : }
1108 :
1109 0 : if (prtlev()>4) { // grad, hess
1110 0 : cout << " grad= " << grad() << endl;
1111 0 : cout << " hess= " << hess() << endl;
1112 : }
1113 :
1114 : } // ibuf
1115 0 : }
1116 :
1117 0 : void VisCalSolver::revert() {
1118 :
1119 0 : if (prtlev()>2) cout << " VCS::revert()" << endl;
1120 :
1121 : // Recall the last decent parameter set
1122 : // TBD: the OK flag?
1123 0 : par()=lastCalPar();
1124 0 : if (svc().solvePol())
1125 0 : srcPar()=lastSrcPar();
1126 0 : }
1127 :
1128 0 : void VisCalSolver::solveGradHess() {
1129 :
1130 0 : if (prtlev()>2) cout << " VCS::solveGradHess()" << endl;
1131 :
1132 : // TBD: explicit option to avoid lmfact?
1133 : // TBD: pointer (or MaskedArray?) optimization?
1134 :
1135 0 : Double lmfact(1.0+lambda());
1136 :
1137 0 : lmfact=2.0;
1138 :
1139 0 : dpar()=Complex(0.0);
1140 0 : for (Int ipar=0; ipar<nCalPar(); ipar++) {
1141 0 : if ( parOK()(ipar) && hess()(ipar)!=0.0) {
1142 : // good hess for this par:
1143 0 : dpar()(ipar) = grad()(ipar)/hess()(ipar);
1144 0 : dpar()(ipar)/=lmfact;
1145 : }
1146 : else {
1147 0 : dpar()(ipar)=0.0;
1148 0 : parOK()(ipar)=false;
1149 : }
1150 : }
1151 :
1152 : // Turn off source pol update for the moment
1153 0 : if (svc().solvePol() && false)
1154 0 : for (Int ipar=nCalPar(); ipar<nTotalPar(); ipar++) {
1155 0 : if ( hess()(ipar)!=0.0) {
1156 : // good hess for this par:
1157 0 : dpar()(ipar) = grad()(ipar)/hess()(ipar);
1158 0 : dpar()(ipar)/=lmfact;
1159 : }
1160 : else {
1161 0 : dpar()(ipar)=0.0;
1162 : }
1163 : }
1164 :
1165 : // Negate (so updatePar() can _add_)
1166 0 : dpar()*=Complex(-1.0f);
1167 :
1168 0 : }
1169 :
1170 0 : void VisCalSolver::updatePar() {
1171 :
1172 0 : if (prtlev()>2) cout << " VCS::updatePar()" << endl;
1173 :
1174 : // if (prtlev()>4) cout << " old =" << par() << endl;
1175 :
1176 : // if (prtlev()>4) cout << " dpar=" << dpar() << endl;
1177 :
1178 :
1179 : // cout << "dCalPar() = " << dCalPar() << endl;
1180 : // cout << "dSrcPar() = " << dSrcPar() << endl;
1181 :
1182 : // Tell svc to update the par
1183 : // (permits svc() to condition the current solutions)
1184 0 : svc().updatePar(dCalPar(),dSrcPar());
1185 :
1186 0 : if (prtlev()>4)
1187 0 : cout << " new =" << endl
1188 0 : << "amp = " << amplitude(par()) << endl
1189 0 : << "pha = " << phase(par()) << endl;
1190 :
1191 0 : }
1192 :
1193 :
1194 0 : void VisCalSolver::optStepSize() {
1195 :
1196 0 : Vector<Double> x2(3,-999.0);
1197 0 : Float step(1.0);
1198 :
1199 : // Starting point is curr chiSq
1200 0 : x2(0)=chiSq();
1201 :
1202 : // take nominal step
1203 0 : par()+=dCalPar();
1204 0 : if (svc().solvePol()) srcPar()+=dSrcPar();
1205 0 : residualate();
1206 0 : chiSquare();
1207 0 : x2(1)=chiSq();
1208 :
1209 : // If nominal step is an improvement...
1210 0 : if (x2(1)<x2(0)) {
1211 :
1212 : // ...double step size until x2 starts increasing
1213 0 : par()=dCalPar(); par()*=Complex(2.0*step); par()+=lastCalPar();
1214 0 : if (svc().solvePol()) { srcPar()=dSrcPar(); srcPar()*=Complex(2.0*step); srcPar()+=lastSrcPar(); };
1215 0 : residualate();
1216 0 : chiSquare();
1217 0 : x2(2)=chiSq();
1218 : // cout << " down: " << step << " " << x2-x2(0) << LogicalArray(x2>=x2(0)) <<endl;
1219 0 : while (x2(2)<x2(1)) { // && step<4.0) {
1220 0 : step*=2.0;
1221 0 : par()=dCalPar(); par()*=Complex(2.0*step); par()+=lastCalPar();
1222 0 : if (svc().solvePol()) { srcPar()=dSrcPar(); srcPar()*=Complex(2.0*step); srcPar()+=lastSrcPar(); };
1223 0 : residualate();
1224 0 : chiSquare();
1225 0 : x2(1)=x2(2);
1226 0 : x2(2)=chiSq();
1227 : // cout << " stretch: " << step << " " << x2-x2(0) << LogicalArray(x2>=x2(0)) <<endl;
1228 :
1229 : }
1230 : }
1231 : // else nominal step too big, so...
1232 : else {
1233 :
1234 : // ... contract by halves until we bracket a minimum
1235 0 : step*=0.5;
1236 0 : par()=dCalPar(); par()*=Complex(step); par()+=lastCalPar();
1237 0 : if (svc().solvePol()) {srcPar()=dSrcPar(); srcPar()*=Complex(step); srcPar()+=lastSrcPar();};
1238 0 : residualate();
1239 0 : chiSquare();
1240 0 : x2(2)=x2(1);
1241 0 : x2(1)=chiSq();
1242 : // cout << " up: " << step << " " << x2-x2(0) << LogicalArray(x2>=x2(0)) <<endl;
1243 0 : while (x2(1)>x2(0)) { // && step>0.125) {
1244 0 : step*=0.5;
1245 0 : par()=dCalPar(); par()*=Complex(step); par()+=lastCalPar();
1246 0 : if (svc().solvePol()) {srcPar()=dSrcPar(); srcPar()*=Complex(step); srcPar()+=lastSrcPar();};
1247 0 : residualate();
1248 0 : chiSquare();
1249 0 : x2(2)=x2(1);
1250 0 : x2(1)=chiSq();
1251 : // cout << " contract: " << step << " " << x2-x2(0) << LogicalArray(x2>=x2(0)) <<endl;
1252 : }
1253 :
1254 : }
1255 :
1256 : // At this point x2(0) > x2(1) < x2(2), so
1257 : // calculate (quadratic) step optimization factor
1258 0 : Double optfactor(0.0);
1259 0 : Double optn(x2(2)-x2(1));
1260 0 : Double optd(x2(0)-2*x2(1)+x2(2));
1261 :
1262 0 : if (abs(optd)>0.0)
1263 0 : optfactor=Double(step)*(1.5-optn/optd);
1264 :
1265 : /*
1266 : cout << "Optimization: "
1267 : << step << " "
1268 : << optfactor << " "
1269 : << x2 << " "
1270 : << "(" << min(amplitude(lastPar())) << ") "
1271 : << max(amplitude(dpar())/amplitude(lastPar()))*180.0/C::pi << " ";
1272 : */
1273 :
1274 0 : par()=lastCalPar();
1275 0 : srcPar()=lastSrcPar();
1276 :
1277 : // Adjust step by the optfactor
1278 0 : if (optfactor>0.0)
1279 0 : dpar()*=Complex(optfactor);
1280 :
1281 : /*
1282 : cout << max(amplitude(dpar())/amplitude(lastPar()))*180.0/C::pi
1283 : << endl;
1284 : */
1285 0 : }
1286 :
1287 0 : void VisCalSolver::optStepSize2() {
1288 :
1289 0 : Vector<Double> x2(3,-999.0);
1290 0 : Float step(1.0);
1291 :
1292 : // Starting point is curr chiSq
1293 0 : x2(0)=chiSq();
1294 :
1295 : // take nominal step
1296 0 : par()+=dCalPar();
1297 0 : if (svc().solvePol()) srcPar()+=dSrcPar();
1298 0 : residualate2();
1299 0 : chiSquare2();
1300 0 : x2(1)=chiSq();
1301 :
1302 : // If nominal step is an improvement...
1303 0 : if (x2(1)<x2(0)) {
1304 :
1305 : // ...double step size until x2 starts increasing
1306 0 : par()=dCalPar(); par()*=Complex(2.0*step); par()+=lastCalPar();
1307 0 : if (svc().solvePol()) { srcPar()=dSrcPar(); srcPar()*=Complex(2.0*step); srcPar()+=lastSrcPar(); };
1308 0 : residualate2();
1309 0 : chiSquare2();
1310 0 : x2(2)=chiSq();
1311 : // cout << " down: " << step << " " << x2-x2(0) << LogicalArray(x2>=x2(0)) <<endl;
1312 0 : while (x2(2)<x2(1)) { // && step<4.0) {
1313 0 : step*=2.0;
1314 0 : par()=dCalPar(); par()*=Complex(2.0*step); par()+=lastCalPar();
1315 0 : if (svc().solvePol()) { srcPar()=dSrcPar(); srcPar()*=Complex(2.0*step); srcPar()+=lastSrcPar(); };
1316 0 : residualate2();
1317 0 : chiSquare2();
1318 0 : x2(1)=x2(2);
1319 0 : x2(2)=chiSq();
1320 : // cout << " stretch: " << step << " " << x2-x2(0) << LogicalArray(x2>=x2(0)) <<endl;
1321 :
1322 : }
1323 : }
1324 : // else nominal step too big, so...
1325 : else {
1326 :
1327 : // ... contract by halves until we bracket a minimum
1328 0 : step*=0.5;
1329 0 : par()=dCalPar(); par()*=Complex(step); par()+=lastCalPar();
1330 0 : if (svc().solvePol()) {srcPar()=dSrcPar(); srcPar()*=Complex(step); srcPar()+=lastSrcPar();};
1331 0 : residualate2();
1332 0 : chiSquare2();
1333 0 : x2(2)=x2(1);
1334 0 : x2(1)=chiSq();
1335 : // cout << " up: " << step << " " << x2-x2(0) << LogicalArray(x2>=x2(0)) <<endl;
1336 0 : while (x2(1)>x2(0)) { // && step>0.125) {
1337 0 : step*=0.5;
1338 0 : par()=dCalPar(); par()*=Complex(step); par()+=lastCalPar();
1339 0 : if (svc().solvePol()) {srcPar()=dSrcPar(); srcPar()*=Complex(step); srcPar()+=lastSrcPar();};
1340 0 : residualate2();
1341 0 : chiSquare2();
1342 0 : x2(2)=x2(1);
1343 0 : x2(1)=chiSq();
1344 : // cout << " contract: " << step << " " << x2-x2(0) << LogicalArray(x2>=x2(0)) <<endl;
1345 : }
1346 :
1347 : }
1348 :
1349 : // At this point x2(0) > x2(1) < x2(2), so
1350 : // calculate (quadratic) step optimization factor
1351 0 : Double optfactor(0.0);
1352 0 : Double optn(x2(2)-x2(1));
1353 0 : Double optd(x2(0)-2*x2(1)+x2(2));
1354 :
1355 0 : if (abs(optd)>0.0)
1356 0 : optfactor=Double(step)*(1.5-optn/optd);
1357 :
1358 : /*
1359 : cout << "Optimization: "
1360 : << step << " "
1361 : << optfactor << " "
1362 : << x2 << " "
1363 : << "(" << min(amplitude(lastPar())) << ") "
1364 : << max(amplitude(dpar())/amplitude(lastPar()))*180.0/C::pi << " ";
1365 : */
1366 :
1367 0 : par()=lastCalPar();
1368 0 : srcPar()=lastSrcPar();
1369 :
1370 : // Adjust step by the optfactor
1371 0 : if (optfactor>0.0)
1372 0 : dpar()*=Complex(optfactor);
1373 :
1374 : /*
1375 : cout << max(amplitude(dpar())/amplitude(lastPar()))*180.0/C::pi
1376 : << endl;
1377 : */
1378 0 : }
1379 :
1380 0 : void VisCalSolver::getErrors() {
1381 :
1382 : // Number of *REAL* dof
1383 : // Int nDOF=2*(nWt()-ntrue(parOK())); // !!!! this is zero for 3 antennas!
1384 0 : Int nDOF=max(2*(nWt()-ntrue(parOK())), 1u);
1385 :
1386 0 : Double k2=chiSq()/Double(nDOF);
1387 :
1388 0 : parErr()=0.0;
1389 : // Vector<Double> snr(nTotalPar(),0.0);
1390 0 : for (Int i=0;i<nCalPar();++i)
1391 0 : if (hess()(i)>0.0) {
1392 0 : parErr()(i)=1.0/sqrt(hess()(i)/k2/2.0); // 2 is from def of Hess!
1393 : // snr(i)=Double(abs(par()(i)))/errs(i);
1394 : }
1395 :
1396 :
1397 : if (false) {
1398 :
1399 : cout << "ChiSq = " << chiSq() << endl;
1400 : cout << "ChiSqV = " << chiSqV() << endl;
1401 : cout << "sumWt = " << sumWt() << endl;
1402 : cout << "nWt = " << nWt()
1403 : << "; nTotalPar() = " << nTotalPar()
1404 : << "; nParOK = " << ntrue(parOK())
1405 : << "; nDOF = " << nDOF
1406 : << endl;
1407 :
1408 : cout << "rChiSq = " << k2 << endl;
1409 : cout << "max(dpar) = " << max(amplitude(dpar())) << endl;
1410 : cout << "Amplitudes = " << amplitude(par()) << endl;
1411 : cout << "Errors = " << mean(parErr()(parOK())) << endl;
1412 :
1413 : }
1414 0 : }
1415 :
1416 :
1417 : } //# NAMESPACE CASA - END
1418 :
|
