Line data Source code
1 : //# UVMod.cc: Implementation of UV Modelling within calibrater
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: casa-feedback@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 : //# $Id$
27 :
28 : #include <synthesis/MeasurementComponents/UVMod.h>
29 :
30 : #include <msvis/MSVis/VisibilityIterator.h>
31 : #include <msvis/MSVis/VisSet.h>
32 : #include <msvis/MSVis/VisBuffer.h>
33 : #include <casacore/casa/Quanta/Quantum.h>
34 : #include <casacore/casa/Quanta/QuantumHolder.h>
35 : #include <casacore/measures/Measures.h>
36 : #include <casacore/measures/Measures/MDirection.h>
37 : #include <casacore/measures/Measures/MEpoch.h>
38 :
39 : #include <casacore/casa/Arrays/ArrayMath.h>
40 : #include <casacore/casa/Arrays/ArrayLogical.h>
41 : #include <casacore/casa/Arrays/ArrayIter.h>
42 : #include <casacore/scimath/Mathematics/MatrixMathLA.h>
43 : #include <casacore/casa/BasicSL/String.h>
44 : #include <casacore/casa/BasicMath/Math.h>
45 : #include <casacore/casa/Utilities/Assert.h>
46 : #include <casacore/casa/Quanta/MVTime.h>
47 : #include <casacore/casa/Exceptions/Error.h>
48 : #include <casacore/casa/OS/Memory.h>
49 : #include <casacore/casa/OS/Path.h>
50 :
51 : #include <components/ComponentModels/ComponentType.h>
52 : #include <components/ComponentModels/Flux.h>
53 : #include <components/ComponentModels/ComponentShape.h>
54 : #include <components/ComponentModels/TwoSidedShape.h>
55 : #include <components/ComponentModels/PointShape.h>
56 : #include <components/ComponentModels/GaussianShape.h>
57 : #include <components/ComponentModels/DiskShape.h>
58 : #include <components/ComponentModels/SpectralModel.h>
59 : #include <components/ComponentModels/ConstantSpectrum.h>
60 : #include <components/ComponentModels/SpectralIndex.h>
61 : #include <components/ComponentModels/SkyCompBase.h>
62 : #include <components/ComponentModels/SkyCompRep.h>
63 : #include <components/ComponentModels/SkyComponent.h>
64 : #include <components/ComponentModels/ComponentList.h>
65 :
66 : #include <sstream>
67 :
68 : #include <casacore/casa/Logging/LogMessage.h>
69 : #include <casacore/casa/Logging/LogSink.h>
70 :
71 :
72 : using namespace casacore;
73 : namespace casa { //# NAMESPACE CASA - BEGIN
74 :
75 :
76 : // **********************************************************
77 : // UVMod Implementations
78 : //
79 :
80 0 : UVMod::UVMod(VisSet& vs) :
81 0 : vs_(NULL),
82 0 : cl_(NULL),
83 0 : svb_(NULL),
84 0 : fitfld_(-1),
85 0 : pc_(),
86 0 : nPar_(0),
87 0 : R_(),
88 0 : dR_(),
89 0 : chiSq_(0.0),
90 0 : lastChiSq_(0.0),
91 0 : sumWt_(0.0),
92 0 : nWt_(0),
93 0 : polWt_(4,false),
94 0 : par_(), lastPar_(),
95 0 : lamb_(0.001), grad_(), lastGrad_(),
96 0 : hess_(), lastHess_(),
97 0 : dpar_(),
98 0 : vary_(),
99 0 : nVary_(0)
100 : {
101 : // cout << "UVMod ctor " << endl;
102 :
103 : // Local copy of VisSet, with correct chunking
104 0 : Block<Int> columns;
105 0 : columns.resize(4);
106 0 : columns[0]=MS::ARRAY_ID;
107 0 : columns[1]=MS::FIELD_ID;
108 0 : columns[2]=MS::DATA_DESC_ID;
109 0 : columns[3]=MS::TIME;
110 :
111 0 : vs_ = new VisSet(vs,columns,3600.0);
112 :
113 : // The VisIter/VisBuffer
114 0 : VisIter& vi(vs_->iter());
115 0 : VisBuffer vb(vi);
116 0 : vi.originChunks();
117 0 : vi.origin();
118 :
119 : // Get phase center & direction ref:
120 0 : pc_=vb.phaseCenter();
121 :
122 : // Check for single field (for now)
123 0 : fitfld()=vb.fieldId();
124 0 : for (vi.originChunks();
125 0 : vi.moreChunks();
126 0 : vi.nextChunk()) {
127 0 : vi.origin();
128 0 : if (vb.fieldId()!=fitfld())
129 0 : throw(AipsError("More than one field Id in selected data"));
130 : }
131 :
132 0 : }
133 :
134 0 : UVMod::~UVMod() {
135 :
136 : // We are only responsible for the ComponentList and VisSet
137 0 : if (cl_) delete cl_; cl_=NULL;
138 0 : if (vs_) delete vs_; vs_=NULL;
139 :
140 0 : }
141 0 : void UVMod::setModel(const ComponentType::Shape type,
142 : const Vector<Double> inpar,
143 : const Vector<Bool> invary) {
144 :
145 : // cout << "UVM::setModel" << endl;
146 : // cout << " type = " << type << endl;
147 : // cout << " inpar = " << inpar << endl;
148 :
149 : // Create empty componentlist
150 0 : cl_ = new ComponentList();
151 :
152 0 : par().resize(inpar.shape());
153 0 : par()=inpar;
154 :
155 : // Make comp
156 0 : SkyComponent comp(type);
157 :
158 0 : switch (type) {
159 0 : case ComponentType::POINT:
160 : {
161 : // A point has 3 pars
162 0 : nPar()=3;
163 :
164 : // Insist par().length()=3
165 0 : if (inpar.nelements()!=3)
166 0 : throw(AipsError("Wrong number of parameters for Point; need 3."));
167 :
168 0 : break;
169 : };
170 :
171 0 : case ComponentType::GAUSSIAN:
172 : case ComponentType::DISK:
173 : {
174 :
175 0 : nPar()=6;
176 :
177 : // Insist par.length()=6
178 0 : if (inpar.nelements()!=6)
179 0 : throw(AipsError("Wrong number of parameters for resolved component; need 6."));
180 :
181 : // Convert pa to radians
182 0 : par()(5)*=(C::pi/180.0);
183 :
184 : // Set shape pars (all in rad)
185 0 : Vector<Double> gpar(3,0.0); // a, b, pa
186 0 : gpar(0)=par()(3)*(C::pi/180.0/60.0/60.0);
187 0 : gpar(1)=gpar(0)*par()(4);
188 0 : gpar(2)=par()(5);
189 0 : comp.shape().setParameters(gpar);
190 :
191 0 : break;
192 0 : };
193 0 : default:
194 : {
195 0 : throw(AipsError("Unrecognized component type in UVMod::setModel."));
196 : break;
197 : }
198 :
199 : }
200 :
201 : // Set I
202 0 : comp.flux().setValue(par()(0));
203 :
204 :
205 : // Add to list
206 0 : cl().add(comp);
207 :
208 :
209 : // Convert direction from arcsec to radians
210 : // par()(1)*=(C::pi/180.0/60.0/60.0);
211 : // par()(2)*=(C::pi/180.0/60.0/60.0);
212 :
213 : // Render direction as an MDirection; use phase-center's DirRef
214 0 : MVDirection mvdir(par()(1)*(C::pi/180.0/60.0/60.0),par()(2)*(C::pi/180.0/60.0/60.0));
215 0 : MDirection dir(mvdir,pc().getRef());
216 :
217 : // Set direction in component
218 0 : skycomp(0).shape().setRefDirection(dir);
219 :
220 : // Handle invary flags:
221 0 : vary().resize(par().shape());
222 :
223 : // Assume varying everything
224 0 : vary()=true;
225 0 : nVary()=nPar();
226 :
227 : // If invary specified, set vary() accordingly
228 0 : for (uInt i=0; i<invary.nelements();i++)
229 0 : if (!invary(i)) {
230 0 : vary()(i) = false;
231 0 : nVary()-=1;
232 : }
233 0 : }
234 :
235 0 : Bool UVMod::modelfit(const Int& maxiter, const String file) {
236 :
237 : // cout << "UVMod::modelfit" << endl;
238 :
239 : // cout << " maxiter = " << maxiter << endl;
240 : // cout << " file = " << file << endl;
241 :
242 0 : LogSink logsink;
243 :
244 : // Local ref to VisIter/VisBuffer
245 0 : VisIter& vi(vs().iter());
246 0 : VisBuffer vb(vi);
247 0 : svb_=&vb;
248 :
249 : // Initialize grad, hess, etc.
250 0 : initSolve();
251 :
252 : {
253 0 : LogMessage message(LogOrigin("UVMod", "solve"));
254 0 : ostringstream o; o<<"Solving for single-component "
255 0 : << skycomp(0).shape().ident();
256 0 : message.message(o);
257 0 : logsink.post(message);
258 0 : }
259 :
260 0 : Int iter(0);
261 : // Begin solution iterations
262 : // Guarantee first pass, which gets initial chi2,
263 : // then evaluate convergence by chi2 comparisons
264 :
265 0 : Bool parOK(true);
266 0 : for (Int validiter=0;validiter<=maxiter;iter++,validiter++) {
267 :
268 : // cout << "iter =" << iter << endl;
269 :
270 : // If we have accumulated Hess/Grad to solve, do so
271 0 : if (iter>0) {
272 :
273 : // If we have done at least one real solve, invoke LM
274 0 : if (iter > 1) {
275 0 : if (parOK && (chiSq()-lastChiSq()) < DBL_EPSILON) {
276 : // if last update succeeded
277 0 : lamb()*=0.5;
278 : // cout << " Good iteration, decreasing lambda: " << lamb() << " " << validiter << endl;
279 0 : lastChiSq()=chiSq();
280 0 : lastPar()=par(); // remember these par,hess,grad,
281 0 : lastGrad()=grad(); // in case we need to resolve
282 0 : lastHess()=hess(); // with new lamb
283 : } else {
284 : // last update failed, so redo with new lamb
285 : // validiter--; // Don't count this iter as valid
286 0 : lamb()*=10.0;
287 :
288 : /*
289 : cout << " Bad iteration, increasing lambda: " << lamb();
290 : if (!parOK) cout << " (bad par update)";
291 : else cout << " (chi2 increased)";
292 : cout << " " << validiter;
293 : cout << endl;
294 : */
295 0 : par()=lastPar(); // recover previous par,hess,grad
296 0 : grad()=lastGrad();
297 0 : hess()=lastHess();
298 : }
299 : } else {
300 : // Remember first par/grad/hess
301 0 : lastChiSq()=chiSq();
302 0 : lastPar()=par(); // remember these par,hess,grad,
303 0 : lastGrad()=grad(); // in case we need to resolve
304 0 : lastHess()=hess(); // with new lamb
305 : }
306 :
307 : // Do the solve
308 0 : solveGradHess();
309 0 : parOK=updPar();
310 : }
311 :
312 : // If pars are ok, calc residuals, etc.
313 0 : if (parOK) {
314 : // zero current chiSq()
315 0 : chiSq()=0.0;
316 0 : sumWt()=0.0;
317 0 : hess()=0.0;
318 0 : grad()=0.0;
319 0 : nWt()=0;
320 :
321 : // Loop over data to accumulate Chi2, Grad, and Hess
322 0 : for (vi.originChunks();
323 0 : vi.moreChunks();
324 0 : vi.nextChunk()) {
325 :
326 : // This version does not pre-apply any calibration!!
327 :
328 : // Loop over VBs:
329 0 : for (vi.origin();
330 0 : vi.more();
331 0 : vi++) {
332 :
333 : // Accumulate chi2 calcuation
334 : // (also accumulates residuals and differentiated residuals)
335 0 : chiSquare();
336 :
337 :
338 : // Accumulate gradients and hessian
339 : // (uses resid and diff'd resid from above)
340 0 : accGradHess();
341 :
342 : }
343 : } // chunks
344 :
345 :
346 :
347 :
348 : // If sum of weights is positive, we can continue with solve
349 0 : if (sumWt()>0) {
350 :
351 : // chiSq()/=sumWt();
352 : // cout << "chiSq(), sumWt() = " << chiSq()<< " " << sumWt() << endl;
353 :
354 : // Report chi2 and pars if a good step
355 0 : if ( (chiSq()-lastChiSq()) < DBL_EPSILON || iter==0) {
356 0 : printPar(validiter);
357 : }
358 : else {
359 0 : validiter--;
360 : // printPar(validiter);
361 : }
362 : } else {
363 : // Escape and report no data!
364 0 : throw(AipsError("Found no data to fit!"));
365 : }
366 :
367 : // Don't get stuck
368 0 : if (iter>100) {
369 0 : cout << "Exceeded maximum iteration limit!" << endl;
370 0 : break;
371 : }
372 :
373 : }
374 : else
375 0 : validiter--;
376 :
377 : }
378 :
379 : //
380 : Double A;
381 0 : A = chiSq()/Double(2*nWt()-nVary());
382 : // grad()/=A;
383 : // hess()/=A;
384 :
385 0 : lamb()=0.0;
386 0 : solveGradHess(true);
387 :
388 : if (false) {
389 : cout << "chiSq() = " << chiSq() << endl;
390 : cout << "sumWt() = " << sumWt() << endl;
391 : cout << "nWt() = " << nWt() << endl;
392 :
393 : cout << "sumWt()/N = " << sumWt()/nWt() << endl;
394 : cout << "<chiSq()> = " << chiSq()/(sumWt()) << endl;
395 :
396 : cout << "A = " << A << endl;
397 :
398 : cout << "grad() = " << grad() << endl;
399 : // cout << "hess() = " << hess() << endl;
400 :
401 : if (nPar()>3) {
402 : dpar()(5)*=(180.0/C::pi);
403 : }
404 :
405 : cout << "dpar() = " << dpar() << endl;
406 :
407 : cout << "cov() = " << hess() << endl;
408 :
409 : Matrix<Double> corr;
410 :
411 : corr = hess();
412 :
413 : for (Int i=0; i<nPar(); i++) {
414 : for (Int j=i; j<nPar(); j++) {
415 : corr(i,j)/=sqrt(hess()(i,i)*hess()(j,j));
416 : }
417 : }
418 :
419 : cout << "corr = " << corr << endl;
420 : }
421 :
422 0 : cout << "If data weights are arbitrarily scaled, the following formal errors" << endl
423 0 : << " will be underestimated by at least a factor sqrt(reduced chi2). If " << endl
424 0 : << " the fit is systematically poor, the errors are much worse." << endl;
425 :
426 0 : Vector<Double> err(nPar(),0.0);
427 0 : for (Int ipar=0; ipar<nPar(); ipar++)
428 0 : if (vary()(ipar))
429 0 : err(ipar)=sqrt(hess()(ipar,ipar));
430 :
431 :
432 0 : cout << "I = " << par()(0) << " +/- " << err(0) << endl;
433 0 : cout << "x = " << par()(1) << " +/- " << err(1) << " arcsec" << endl;
434 0 : cout << "y = " << par()(2) << " +/- " << err(2) << " arcsec" << endl;
435 0 : if (nPar()>3) {
436 0 : cout << "a = " << par()(3) << " +/- " << err(3) << " arcsec" << endl;
437 0 : cout << "r = " << par()(4) << " +/- " << err(4) << endl;
438 0 : cout << "p = " << par()(5)*180.0/C::pi << " +/- " << err(5)*180.0/C::pi << " deg" << endl;
439 : }
440 :
441 : // } // (false)
442 :
443 : // Shift pa back to deg
444 0 : if (par().nelements()>3) {
445 0 : par()(5)*=(180.0/C::pi);
446 : }
447 :
448 : // Shift model to phase center of selected field
449 0 : MDirection newdir;
450 0 : newdir=pc();
451 0 : newdir.shift(skycomp(0).shape().refDirection().getValue(),true);
452 0 : skycomp(0).shape().setRefDirection(newdir);
453 :
454 : // Export componentlist to file, if specified
455 0 : if (file.length()>0) {
456 0 : Path path(file);
457 0 : cout << "Writing componentlist to file: "
458 0 : << path.absoluteName()
459 0 : << endl;
460 0 : cl().rename(path);
461 0 : }
462 :
463 0 : return true;
464 :
465 0 : }
466 :
467 0 : void UVMod::initSolve() {
468 :
469 : // cout << "UVM::initSolve" << endl;
470 :
471 : // Chi2 and weights
472 0 : chiSq()=0.0;
473 0 : sumWt()=0.0;
474 0 : nWt()=0;
475 :
476 0 : lastPar().resize(nPar());
477 :
478 : // Gradient and Hessian
479 0 : grad().resize(nPar());
480 0 : grad()=0.0;
481 :
482 0 : lastGrad().resize(nPar());
483 0 : lastGrad()=0.0;
484 :
485 0 : hess().resize(nPar(),nPar());
486 0 : hess()=0.0;
487 :
488 0 : lastHess().resize(nPar(),nPar());
489 0 : lastHess()=0.0;
490 :
491 0 : dpar().resize(nPar());
492 0 : dpar()=0.0;
493 :
494 0 : lamb()=0.001;
495 :
496 0 : }
497 :
498 0 : void UVMod::chiSquare() {
499 :
500 : // cout << "UVM::chiSquare" << endl;
501 :
502 : // Get residuals w.r.t. current parameters
503 0 : residual();
504 :
505 : // Pointer access to vb elements
506 0 : const Bool* flr=&svb().flagRow()(0);
507 0 : const Bool* fl= &svb().flag()(0,0);
508 0 : const Int* a1= &svb().antenna1()(0);
509 0 : const Int* a2= &svb().antenna2()(0);
510 0 : const Float* wt= &svb().weight()(0);
511 :
512 :
513 0 : Int nCorr = R().shape()(0);
514 :
515 : DComplex* res;
516 :
517 : // Loop over row/channel
518 0 : Int irow(0),ich(0),icorr(0);
519 0 : for (irow=0;irow<svb().nRow();irow++,flr++,a1++,a2++,wt++) {
520 0 : if (!(*flr) && (*a1)!=(*a2)) { // not flagrow nor AC
521 0 : fl=&svb().flag()(0,irow);
522 0 : for (ich=0;ich<svb().nChannel();ich++,fl++) {
523 0 : if (!(*fl)) {
524 :
525 0 : res=&R()(0,ich,irow);
526 0 : for (icorr=0; icorr<nCorr; icorr++,res++) {
527 :
528 0 : if ( polWt()(icorr) ) {
529 0 : chiSq() += Double(*wt)*real((*res)*conj(*res));
530 0 : sumWt() += Double(*wt);
531 0 : nWt() += 1;
532 :
533 : /*
534 : cout << irow << " " << *a1 << " " << *a2 << " ";
535 : cout << "*res = " << *res << " ";
536 : cout << innerProduct((*res),(*res)) << " ";
537 : cout << chiSq()/sumWt();
538 : cout << endl;
539 : */
540 :
541 : // chiSq() += real(innerProduct((*res),(*res)));
542 : // sumWt() += 1.0;
543 : }
544 : }
545 :
546 : }
547 : }
548 : }
549 : }
550 :
551 : // cout << "End of UVM::chiSquare" << chiSq() << endl;
552 :
553 0 : }
554 :
555 :
556 : // Vobs - (M)*Vmod, - (dM)*Vmod
557 0 : void UVMod::residual() {
558 :
559 : // cout << "UVM::residual" << endl;
560 :
561 : // Shapes
562 0 : Int nRow= svb().nRow();
563 0 : Int nChan=svb().nChannel();
564 :
565 : // Data Access
566 0 : const Bool* flagR= &svb().flagRow()(0);
567 0 : const Bool* flag= &svb().flag()(0,0);
568 0 : const Int* a1= &svb().antenna1()(0);
569 0 : const Int* a2= &svb().antenna2()(0);
570 0 : const RigidVector<Double,3>* uvw = &svb().uvw()(0);
571 :
572 : // Prepare residuals array
573 0 : Int nCorr = svb().correctedVisCube().shape()(0);
574 0 : R().resize(nCorr,nChan,nRow);
575 0 : convertArray(R(),svb().correctedVisCube());
576 :
577 : // Zero cross-hands...
578 : // TBD
579 :
580 : // Prepare differentials array
581 0 : dR().resize(IPosition(4,nCorr,nPar(),nChan,nRow));
582 0 : dR()=DComplex(0.0);
583 :
584 : // Calculate wave number per frequency
585 0 : Vector<Double> freq = svb().frequency();
586 0 : Vector<Double> waveNum = C::_2pi*freq/C::c;
587 :
588 : // cout << "waveNum = " << waveNum << endl;
589 :
590 : // cout << "Polarizations = " << svb().corrType() << endl;
591 :
592 :
593 0 : Vector<Int> corridx = svb().corrType();
594 :
595 : // Ensure component is in correct pol
596 0 : ComponentType::Polarisation pol(ComponentType::CIRCULAR);
597 0 : if (svb().polFrame()==0) {
598 0 : pol=ComponentType::CIRCULAR;
599 0 : corridx-=5;
600 : }
601 0 : else if (svb().polFrame()==1) {
602 0 : pol=ComponentType::LINEAR;
603 0 : corridx-=9;
604 : }
605 :
606 : // Handle polarization selections (parallel hands only, for now)
607 0 : polWt().resize(nCorr);
608 0 : polWt()=false;
609 0 : polWt() = (corridx==0 || corridx==3);
610 :
611 0 : skycomp(0).flux().convertPol(pol);
612 :
613 : // The per-row model visibility
614 0 : Vector<DComplex> M(4);
615 :
616 : // The comp flux
617 0 : Vector<DComplex> F(4);
618 0 : F = skycomp(0).flux().value();
619 :
620 : // The Direction visibility, derivatives
621 0 : DComplex D, dphdx, dphdy;
622 :
623 : // short-hand access to pars
624 0 : Vector<Double> p;
625 0 : p.reference(par());
626 :
627 : // Utility
628 0 : Vector<DComplex> dMdG;
629 : Double cospa;
630 : Double sinpa;
631 0 : Double dGda(0.0), G(0.0), g1(0.0),g2(0.0);
632 0 : Double pi2a2(0.0);
633 0 : if (nPar()>3)
634 0 : pi2a2=C::pi*C::pi*p(3)*p(3);
635 :
636 : // Pointer access to R(), dR();
637 0 : DComplex *res = &R()(0,0,0);
638 0 : DComplex *dres = &dR()(IPosition(4,0,0,0,0));
639 :
640 : // iterate rows
641 0 : for (Int row=0; row<nRow; row++,flagR++,uvw++,a1++,a2++) {
642 :
643 0 : if (!*flagR && (*a1)!=(*a2)) { // not flagrow nor AC
644 :
645 0 : flag = &svb().flag()(0,row);
646 0 : for (Int chn=0; chn<nChan; chn++,flag++) {
647 :
648 : // if this data channel unflagged
649 0 : if (!*flag) {
650 :
651 0 : Vector<Double> uvw2(uvw->vector());
652 : // The model visbility vector
653 0 : M=skycomp(0).visibility(uvw2,freq(chn)).value();
654 :
655 : // Scale uvw2 to 1/arcsec (dir pars are in arcsec)
656 0 : uvw2/=(648000/C::pi);
657 :
658 : // Direction phase factors (OPTIMIZED?)
659 0 : Double phase=waveNum(chn)*( uvw2(0)*p(1)+uvw2(1)*p(2) );
660 0 : D=DComplex(cos(phase),sin(phase));
661 0 : dphdx=DComplex(0.0,uvw2(0)*waveNum(chn));
662 0 : dphdy=DComplex(0.0,uvw2(1)*waveNum(chn));
663 :
664 : // apply direction phase
665 0 : M*=D;
666 :
667 : // Resolved component geometry derivatives
668 0 : if (nPar()>3) {
669 :
670 0 : cospa=cos(p(5));
671 0 : sinpa=sin(p(5));
672 0 : g1=waveNum(chn)*(uvw2(0)*cospa - uvw2(1)*sinpa)/C::_2pi;
673 0 : g2=waveNum(chn)*(uvw2(0)*sinpa + uvw2(1)*cospa)/C::_2pi;
674 0 : dGda=C::pi*sqrt(g1*g1*p(4)*p(4) + g2*g2);
675 0 : G=p(3)*dGda;
676 :
677 : // Gaussian-specific, for now
678 0 : switch (skycomp(0).shape().type()) {
679 0 : case ComponentType::GAUSSIAN:
680 : {
681 0 : dMdG=M*DComplex(-G/2.0/C::ln2);
682 0 : break;
683 : }
684 0 : case ComponentType::DISK:
685 : {
686 0 : dMdG=M*DComplex(-2.0*jn(2,G)/G);
687 0 : break;
688 : }
689 0 : default:
690 : {
691 0 : break;
692 : }
693 : }
694 : }
695 :
696 : // Fill R, dR by element:
697 0 : res = &R()(0,chn,row);
698 0 : dres = &dR()(IPosition(4,0,0,chn,row));
699 0 : Int cidx(0);
700 0 : for (Int icorr=0; icorr<nCorr; icorr++,res++,dres++) {
701 :
702 0 : if (polWt()(icorr)) {
703 :
704 : // Re-index correlations in model
705 : // (copes with partial or mis-ordered data correlations)
706 0 : cidx=corridx(icorr);
707 :
708 : // The residual itself
709 0 : *(res) -= M(cidx);
710 :
711 : // Flux derivative
712 0 : *(dres) = -M(cidx)/F(cidx);
713 :
714 : // Direction derivatives
715 0 : *(dres+1*nCorr) = -M(cidx)*dphdx;
716 0 : *(dres+2*nCorr) = -M(cidx)*dphdy;
717 :
718 : // If fitting resolved shapes:
719 0 : if (nPar()>3) {
720 :
721 : // 2D-comp-specific par derivatives
722 :
723 : // a
724 0 : *(dres+3*nCorr) = -dMdG(cidx)*dGda;
725 :
726 : // r = b/a
727 0 : *(dres+4*nCorr) = -dMdG(cidx)*(pi2a2*p(4)*g1*g1/G);
728 :
729 : // pa
730 0 : *(dres+5*nCorr) = -dMdG(cidx)*(pi2a2*(1.0-p(4)*p(4))*g1*g2/G);
731 : }
732 : }
733 : }
734 0 : } // !*flag
735 : } // chn
736 : } // !*flagR
737 : } // row
738 :
739 : // cout << "End of UVM::residual" << endl;
740 :
741 :
742 0 : }
743 :
744 0 : void UVMod::accGradHess() {
745 :
746 : // Assumes we've already calculated R() & dR() (via chiSquare/residual)
747 :
748 : // cout << "UVM::accGradHess" << endl;
749 :
750 0 : const Bool* flagR= &svb().flagRow()(0);
751 0 : const Bool* flag= &svb().flag()(0,0);
752 0 : const Int* a1= &svb().antenna1()(0);
753 0 : const Int* a2= &svb().antenna2()(0);
754 0 : const Float* wt= &svb().weight()(0);
755 :
756 0 : Int nCorr = R().shape()(0);
757 :
758 : DComplex *res;
759 : DComplex *dres1, *dres2;
760 :
761 : // grad()=0.0;
762 : // hess()=0.0;
763 :
764 : // cout << endl;
765 :
766 0 : for (Int irow=0;irow<svb().nRow();irow++,a1++,a2++,flagR++,wt++) {
767 :
768 0 : if (!(*flagR)) {
769 :
770 0 : flag= &svb().flag()(0,irow);
771 0 : for (Int ich=0;ich<svb().nChannel();ich++,flag++) {
772 :
773 0 : if (!(*flag)) {
774 :
775 0 : res = &R()(0,ich,irow);
776 0 : for (Int icorr=0; icorr<nCorr; icorr++,res++) {
777 :
778 0 : if ( polWt()(icorr) ) {
779 0 : dres1 = &dR()(IPosition(4,icorr,0,ich,irow));
780 0 : for (Int ipar0=0;ipar0<nPar();ipar0++,dres1+=nCorr) {
781 :
782 : // Only if this parameter varies, calc grad/hess
783 0 : if (vary()(ipar0)) {
784 0 : grad()(ipar0) += (Double(*wt)*real( (*res) * conj(*dres1) ) );
785 0 : hess()(ipar0,ipar0) += (Double(*wt)*real( (*dres1) * conj(*dres1) ) );
786 :
787 0 : dres2 = dres1+nCorr;
788 0 : for (Int ipar1=ipar0+1;ipar1<nPar();ipar1++,dres2+=nCorr) {
789 0 : if (vary()(ipar1))
790 0 : hess()(ipar0,ipar1) += (Double(*wt)*real( (*dres1) * conj(*dres2) ) );
791 : }
792 : }
793 : }
794 : }
795 : }
796 : }
797 :
798 : }
799 : }
800 : }
801 :
802 : // cout << "grad() = " << grad().column(1) << endl;
803 : // cout << "hess() = " << hess().xyPlane(1) << endl;
804 :
805 :
806 :
807 0 : }
808 :
809 0 : void UVMod::solveGradHess(const Bool& doCovar) {
810 :
811 : // cout << "UVM::solveGradHess" << endl;
812 :
813 : // Treat diagonal
814 0 : for (Int ipar=0; ipar<nPar(); ipar++) {
815 : // Ensure non-zero diag
816 0 : if (hess()(ipar,ipar)==0.0) // corresponds to vary()(ipar)=false
817 0 : hess()(ipar,ipar)=1.0;
818 :
819 : // apply lamb() to diag (if covar matrix not requested)
820 0 : if (!doCovar) hess()(ipar,ipar)*=(1.0+lamb());
821 : }
822 :
823 : // cout << "grad() = " << grad() << endl;
824 : // cout << "hess() = " << hess() << endl;
825 :
826 :
827 0 : char uplo('U');
828 0 : Int n(nPar());
829 0 : Int nrhs(1);
830 : Bool deleteaa, deletebb;
831 0 : Int lda(n);
832 0 : Int ldb(n);
833 : Int info;
834 : Double *aa;
835 : Double *bb;
836 :
837 0 : if (sumWt() > 0.0) {
838 :
839 0 : aa = hess().getStorage(deleteaa);
840 0 : bb = grad().getStorage(deletebb);
841 :
842 0 : posv(&uplo, &n, &nrhs, aa, &lda, bb, &ldb, &info);
843 :
844 0 : if (doCovar) potri(&uplo, &n, aa, &lda, &info);
845 :
846 : // if (info!=0) cout << "info = " << info << endl;
847 :
848 0 : hess().putStorage(aa,deleteaa);
849 0 : grad().putStorage(bb,deletebb);
850 : }
851 0 : dpar()=grad();
852 :
853 : // cout << "dpar() = " << dpar() << endl;
854 :
855 0 : }
856 :
857 0 : Bool UVMod::updPar() {
858 :
859 : // cout << "UVM::updPar" << endl;
860 :
861 : // cout << endl << "dpar() = " << dpar() << endl;
862 :
863 :
864 : // Handle "unphysical" updates:
865 :
866 0 : if (nPar() > 3 ) {
867 0 : if (dpar()(3) > par()(3)) // Keep size positive
868 0 : dpar()(3) = 0.9*par()(3);
869 0 : if (dpar()(4) > par()(4)) // Keep axial ratio positive
870 0 : dpar()(4) = 0.9*par()(4);
871 : }
872 :
873 : // Update pars
874 0 : par()-=dpar();
875 :
876 0 : if (nPar()>3 && par()(4) > 1.0) par()(4)=1.0; // Keep axial ration <= 1
877 :
878 :
879 : // Set current pars in the componentlist
880 0 : return setCompPar();
881 :
882 : }
883 :
884 0 : Bool UVMod::setCompPar() {
885 :
886 : // Update I flux:
887 0 : skycomp(0).flux().convertPol(ComponentType::STOKES);
888 0 : skycomp(0).flux().setValue(par()(0));
889 :
890 0 : Double rad2as(180.0*60.0*60.0/C::pi);
891 :
892 : // Update (relative) direction:
893 0 : MVDirection mvdir(par()(1)/rad2as,par()(2)/rad2as);
894 0 : MDirection dir(mvdir);
895 0 : skycomp(0).shape().setRefDirection(dir);
896 :
897 : // Set shape pars
898 0 : if (skycomp(0).shape().nParameters() > 0) {
899 :
900 : // Keep pa within a quarter cycle of zero
901 0 : while (par()(5) > C::pi/2.0) par()(5)-=C::pi;
902 0 : while (par()(5) < -C::pi/2.0) par()(5)+=C::pi;
903 :
904 0 : Vector<Double> gpar(skycomp(0).shape().parameters()); // a, b, pa
905 : // cout << "gpar = " << gpar << endl;
906 0 : gpar(0)=par()(3)/rad2as; // a (rad)
907 0 : gpar(1)=gpar(0)*par()(4); // b = a*r (rad)
908 0 : gpar(2)=par()(5); // pa (rad)
909 :
910 :
911 : // Now set new pars
912 : // cout << "gpar = " << gpar << endl;
913 : try {
914 0 : skycomp(0).shape().setParameters(gpar);
915 0 : } catch (AipsError& x) {
916 0 : cout << " This should never happen now: " << x.getMesg() << endl;
917 0 : return false;
918 0 : }
919 0 : }
920 0 : return true;
921 0 : }
922 :
923 0 : void UVMod::printPar(const Int& iter) {
924 :
925 : // cout << "UVM::printPar" << endl;
926 :
927 : // Ensure component is in Stokes:
928 0 : skycomp(0).flux().convertPol(ComponentType::STOKES);
929 :
930 : // Extract flux:
931 0 : Vector<Double> f;
932 0 : skycomp(0).flux().value(f);
933 :
934 0 : if (iter==0) {
935 0 : cout << "There are "
936 0 : << 2*nWt() << " - "
937 0 : << nVary() << " = "
938 0 : << 2*nWt()-nVary() << " degrees of freedom."
939 0 : << endl;
940 : }
941 :
942 0 : cout << " iter=" << iter << ": ";
943 :
944 : // cout << " lastchi2=" << lastChiSq() << " ";
945 :
946 0 : if (iter>0 && (chiSq()-lastChiSq())>DBL_EPSILON ) cout << "(";
947 0 : cout << " reduced chi2=" << chiSq()/Double(2*nWt()-nVary());
948 0 : if (iter>0 && (chiSq()-lastChiSq())>DBL_EPSILON ) cout << ")";
949 :
950 : // cout << "(" << lamb() << ")";
951 :
952 0 : cout << ": I=" << f(0)
953 0 : << ", dir="
954 0 : << skycomp(0).shape().refDirection().getAngle(Unit("arcsec"));
955 :
956 0 : if (skycomp(0).shape().nParameters()>0) {
957 0 : Vector<Double> gpar(skycomp(0).shape().parameters()); // a, b, pa
958 0 : cout << ", shape=["
959 0 : << gpar(0)*180.0*60.0*60.0/C::pi << ", " // a (arcsec)
960 0 : << gpar(1)/gpar(0) << ", " // r
961 0 : << gpar(2)*180.0/C::pi << "]"; // pa (deg)
962 0 : }
963 0 : cout << endl;
964 :
965 0 : }
966 :
967 : } //# NAMESPACE CASA - END
968 :
|
