Line data Source code
1 : //# ConjugateGradientSolver.cc: Implementation of an iterative ConjugateGradientSolver
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 adressed 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 : //# $Id: ConjugateGradientSolver.cc,v 1.4 2006/05/02 15:59:22 sbhatnag Exp $
28 :
29 : #include <synthesis/MeasurementComponents/SteepestDescentSolver.h>
30 : #include <synthesis/TransformMachines/Utils.h>
31 : #include <casacore/casa/Arrays/MatrixMath.h>
32 : #include <limits>
33 :
34 : using namespace casacore;
35 : namespace casa {
36 :
37 0 : SteepestDescentSolver::SteepestDescentSolver(Int nParams,Vector<Int> polMap,
38 0 : Int nIter, Double tol):
39 0 : Iterate(), polMap_p(polMap)
40 : {
41 0 : setNumberIterations(nIter);
42 0 : setTolerance(tol);
43 0 : setMaxParams(nParams);
44 0 : };
45 : //
46 : //---------------------------------------------------------------------------
47 : // Get a vector of visibilities from all baselines with antenna
48 : // whichAnt. It returns N-1 values (fills in the complex conjugate
49 : // part as well).
50 : //
51 0 : Vector<Complex> SteepestDescentSolver::getVj(const VisBuffer& vb, Int NAnt, Int whichAnt,
52 : Int whichPol,
53 : Double& sumWt,Int negate,Int weighted)
54 : {
55 0 : Vector<Int> ant1, ant2;
56 0 : Vector<Complex> Vj;
57 0 : Vector<Int> nPoints;
58 0 : ant1 = vb.antenna1();
59 0 : ant2 = vb.antenna2();
60 0 : Vj.resize(NAnt);
61 0 : nPoints.resize(NAnt);
62 0 : Int J=0,N;
63 0 : IPosition ndx(3,0);
64 0 : Vj = 0;
65 0 : nPoints = 0;
66 0 : N = vb.nRow();
67 : Double wt;
68 0 : sumWt = 0;
69 : //
70 : // It will be most useful to use
71 : // both RR and LL, and both cases are useful: (1) compute joind RR
72 : // and LL chisq, and (2) compute separate Chisq for RR and LL.
73 : // The CalTable stuff below all this is however not yet capable of
74 : // writing a caltable with 2 real parameters for 2 polarizations.
75 : //
76 : // Currently computing chisq for only the first polarization in
77 : // vis. ArrayColumn from the MS.
78 : //
79 0 : ndx=0;ndx(0)=whichPol;
80 : Int i;
81 : //#pragma omp parallel default(shared) private(i)
82 : {
83 : //#pragma omp for
84 : // for(ndx(2);ndx(2)<N;ndx(2)++)
85 0 : for(i=ndx(2);i<N;i++)
86 : {
87 0 : ndx(2)=i;
88 0 : if ((!vb.flagRow()(ndx(2))) &&
89 0 : (!vb.flag()(0,ndx(2))) &&
90 : // (!vb.flagCube()(ndx)) &&
91 0 : ((ant1[ndx(2)] != ant2[ndx(2)]) &&
92 0 : ((ant1[ndx(2)] == whichAnt) ||
93 0 : (ant2[ndx(2)] == whichAnt))
94 : ))
95 : {
96 0 : wt = vb.weight()(ndx(2));
97 0 : sumWt += wt;
98 :
99 0 : if (!weighted) wt = 1.0;
100 :
101 0 : J=(ant1[ndx[2]]!=whichAnt)?ant1[ndx(2)]:ant2[ndx(2)];
102 0 : nPoints[J]++;
103 0 : if (ant1[ndx(2)] > ant2[ndx(2)])
104 0 : Vj[J] += Complex(wt,0)*vb.modelVisCube()(ndx);
105 : else
106 0 : if (negate)
107 0 : Vj[J] += -Complex(wt,0)*(vb.modelVisCube()(ndx));
108 : else
109 0 : Vj[J] += Complex(wt,0)*conj(vb.modelVisCube()(ndx));
110 0 : J++;
111 : }
112 : }
113 : }
114 : //#pragma omp parallel default(shared) private(i)
115 : //#pragma omp for
116 0 : for(i=0;i<NAnt;i++)
117 0 : if (nPoints[i] <= 0) Vj[i] = 0.0;
118 :
119 0 : return Vj;
120 : }
121 : //
122 : //---------------------------------------------------------------------------
123 : // Compute the penalty function (also called the Goodness-of-fit criteria).
124 : // For us, its the Chi-square function.
125 : //
126 0 : Double SteepestDescentSolver::getGOF(const VisBuffer& residual, Int& whichPol, Double& sumWt, const char* /*msg*/)
127 : {
128 0 : Double Chisq=0.0;
129 0 : Int nRow=residual.nRow();
130 0 : IPosition vbShape(residual.modelVisCube().shape());
131 : //Int nPoints;
132 : //nPoints=0;
133 0 : sumWt=0.0;
134 :
135 0 : IPosition ndx(3,0);
136 0 : ndx(0)=whichPol;
137 0 : for (ndx(2)=0;ndx(2)<nRow;ndx(2)++)
138 : {
139 0 : if (
140 : // (!residual.flagCube()(ndx)) &&
141 0 : (!residual.flag()(0,ndx(2))) &&
142 0 : (!residual.flagRow()(ndx(2))) &&
143 0 : (residual.antenna1()(ndx(2)) != residual.antenna2()(ndx(2)))
144 : )
145 : {
146 0 : DComplex Vis;
147 : Float wt;
148 :
149 0 : Vis = residual.modelVisCube()(ndx);
150 0 : wt = residual.weight()(ndx(2));
151 : /*
152 : cout << ndx(2) << " "
153 : << Vis << " "
154 : << wt << " "
155 : << residual.antenna1()(ndx(2)) << "-" << residual.antenna2()(ndx(2)) << " "
156 : << residual.flagCube()(ndx) << " "
157 : << residual.flag()(0,ndx(2)) << " "
158 : << residual.flagRow()(ndx(2)) << " "
159 : << endl;
160 : */
161 0 : sumWt += Double(wt);
162 0 : Chisq += Double(wt*real(Vis*conj(Vis)));
163 : // nPoints++;
164 : /*
165 : if (isnan(Chisq))
166 : {
167 : ostringstream os;
168 : os << "###Error: Chisq evaluates to a NaN. SumWt = " << sumWt
169 : << ndx(2) << " " << Vis << " " << wt << " "
170 : << residual.antenna1()(ndx(2)) << " " << residual.antenna2()(ndx(2));
171 : cout << residual.modelVisCube() << endl;
172 : throw(AipsError(os.str().c_str()));
173 : }
174 : */
175 : }
176 : }
177 : // exit(0);
178 : // if (nPoints>0) sumWt /= nPoints;
179 : //cout << "CHISQ CHISQ = " << Chisq/sumWt << " " << Chisq << " " << sumWt << endl;
180 : /*
181 : ndx(0)=whichPol;
182 : for (ndx(2)=0;ndx(2)<nRow;ndx(2)++)
183 : // for (ndx(2)=0;ndx(2)<2;ndx(2)++)
184 : {
185 : Int i = ndx(2),p=ndx(0);
186 : if (
187 : (!residual.flagRow()(ndx(2))) &&
188 : (!residual.flag()(0,ndx(2))) &&
189 : // (!residual.flagCube()(ndx)) &&
190 : (residual.antenna1()(ndx(2)) != residual.antenna2()(ndx(2)))
191 : )
192 : cout << "SDS: Residual: " << msg << " " << i << " " << whichPol
193 : << " " << getCurrentTimeStamp(residual)-4.46635e9
194 : << " " << residual.modelVisCube()(0,0,i)
195 : << " " << residual.modelVisCube()(1,0,i)
196 : << " " << residual.visCube()(0,0,i)
197 : << " " << residual.visCube()(1,0,i)
198 : << " " << residual.antenna1()(i)<< "-" << residual.antenna2()(i)
199 : << " " << residual.flag()(0,i)
200 : << " " << residual.flagRow()(i)
201 : << " " << residual.flagCube()(p,0,i)
202 : << " " << residual.weight()(i)
203 : << endl;
204 : }
205 : cout << "Chisq = " << Chisq << " " << sumWt << " " << Chisq/sumWt << endl;
206 : */
207 : // exit(0);
208 0 : if (sumWt > 0) return Chisq/sumWt;
209 0 : else return 0.0;
210 : }
211 : //
212 : //---------------------------------------------------------------------------
213 : //Given the VisEquation, this iteratively sovles for the parameters
214 : //of the supplied VisJones for the time-stamp given by SlotNo. nAnt
215 : //is the number of antennas per time stamp.
216 : //
217 0 : Double SteepestDescentSolver::solve(VisEquation& ve, EPJones& epj, VisBuffer& vb,
218 : Int nAnt, Int SlotNo)
219 : {
220 0 : Cube<Float> oldOffsets;
221 0 : Vector<Complex> ResidualVj, dAzVj, dElVj;
222 : Double Chisq0,Chisq,sumWt;
223 : //Double Time;
224 : //static Double Time0;
225 : Double AzHDiag,ElHDiag;
226 0 : Timer timer;
227 : Double localGain;
228 0 : Int iPol=0;
229 0 : IPosition ndx(2,2,nAnt);
230 :
231 0 : localGain = gain();
232 0 : ResidualVj.resize(nAnt);
233 0 : dAzVj.resize(nAnt);
234 0 : dElVj.resize(nAnt);
235 :
236 0 : ndx(0)=1;
237 0 : Chisq0=0.0;
238 0 : epj.guessPar(vb);
239 0 : epj.solveRPar() = -0.0001;
240 0 : residual_p=vb;
241 0 : residual_p.visCube().resize(vb.visCube().shape());
242 0 : residual_p.modelVisCube().resize(vb.modelVisCube().shape());
243 :
244 0 : residual_p.visCube() = vb.visCube();
245 :
246 0 : ve.diffResiduals(residual_p, gradient0_p, gradient1_p,flags);
247 :
248 : try
249 : {
250 0 : Chisq = getGOF(residual_p,iPol,sumWt);
251 : }
252 0 : catch (AipsError &x)
253 : {
254 0 : cout << x.getMesg() << endl;
255 0 : return -1;
256 : }
257 :
258 0 : if (sumWt == 0) return -Chisq;
259 :
260 : //Time = getCurrentTimeStamp(residual_p);
261 : //if (SlotNo == 0) Time0 = Time;
262 :
263 0 : timer.mark();
264 : Int iter;
265 : //Bool Converged=false;
266 0 : Vector<Bool> antFlagged;
267 0 : antFlagged.resize(nAnt);
268 0 : logIO() << LogOrigin("PointingCal","SDS::solve()")
269 : << "Solution interval = " << SlotNo
270 : << " Initial Chisq = " << Chisq
271 0 : << " SumOfWt = " << sumWt << LogIO::NORMAL3 << LogIO::POST;
272 :
273 0 : for (iter=0;iter<numberIterations();iter++)
274 : {
275 0 : oldOffsets = epj.solveRPar();//Guess;
276 : {
277 0 : for(Int ant=0;ant<nAnt;ant++)
278 : {
279 0 : ResidualVj = getVj(residual_p, nAnt,ant,iPol,sumWt,0,0);//Get a weighted list
280 0 : dAzVj = getVj(gradient0_p,nAnt,ant,iPol,sumWt,0,1);// Get an un-weighted list
281 0 : dElVj = getVj(gradient1_p,nAnt,ant,iPol,sumWt, 0,1);
282 :
283 0 : antFlagged[ant]=false;
284 0 : if (sumWt > 0)
285 : {
286 0 : ndx(1)=ant;
287 :
288 : Double coVar1, coVar2;
289 0 : coVar1 = coVar2 = 0;
290 0 : if (sumWt > 0) coVar1 = real(innerProduct(dAzVj,dAzVj))/(sumWt*sumWt);
291 0 : if (sumWt > 0) coVar2 = real(innerProduct(dElVj,dElVj))/(sumWt*sumWt);
292 : // if (sumWt > 0) coVar1 = real(innerProduct(dAzVj,dAzVj));
293 : // if (sumWt > 0) coVar2 = real(innerProduct(dElVj,dElVj));
294 0 : AzHDiag = ElHDiag = 0;
295 0 : AzHDiag = 1.0/coVar1;
296 0 : ElHDiag = 1.0/coVar2;
297 :
298 0 : epj.solveRPar()(0,iPol,ant) = epj.solveRPar()(0,iPol,ant)-2*AzHDiag*localGain*
299 0 : real(innerProduct(ResidualVj,dAzVj))/sumWt;
300 0 : epj.solveRPar()(1,iPol,ant) = epj.solveRPar()(1,iPol,ant)-2*ElHDiag*localGain*
301 0 : real(innerProduct(ResidualVj,dElVj))/sumWt;
302 :
303 : // epj.solveRPar()(0,0,ant) = epj.solveRPar()(0,0,ant)-2*AzHDiag*localGain*
304 : // real(innerProduct(ResidualVj,dAzVj));
305 : // epj.solveRPar()(1,0,ant) = epj.solveRPar()(1,0,ant)-2*ElHDiag*localGain*
306 : // real(innerProduct(ResidualVj,dElVj));
307 : }
308 : else
309 0 : antFlagged[ant]=true;
310 :
311 : }
312 : //
313 : // Compute the residuals and the gradients with the
314 : // updated solutions and check for convergence.
315 : //
316 : // epj.setAntPar(SlotNo,Guess);
317 0 : ve.diffResiduals(residual_p,gradient0_p,gradient1_p,flags);
318 0 : Chisq0 = Chisq;
319 : // Chisq = getGOF(residual,iPol,sumWt)/sumWt;
320 : try
321 : {
322 0 : Chisq = getGOF(residual_p,iPol,sumWt);
323 : }
324 0 : catch (AipsError &x)
325 : {
326 0 : logIO() << LogOrigin("PointingCal","SDS::solve()")
327 : << x.getMesg()
328 0 : << LogIO::POST;
329 0 : return -1;
330 : }
331 : }
332 : Double dChisq;
333 0 : dChisq = (Chisq0-Chisq);
334 : // cout << iter << " " << Chisq0 << " " << Chisq << " " << dChisq << endl;
335 0 : if ((fabs(dChisq) < tolerance()))
336 : {
337 : //Converged=true;
338 0 : break;
339 : }
340 0 : if ((dChisq < tolerance()))// && (iter > 0))
341 : {
342 0 : localGain /= 2.0;
343 0 : epj.solveRPar() = oldOffsets;
344 0 : ve.diffResiduals(residual_p,gradient0_p,gradient1_p,flags);
345 : // Chisq = getGOF(residual,iPol,sumWt)/sumWt;
346 : try
347 : {
348 0 : Chisq = getGOF(residual_p,iPol,sumWt);
349 : }
350 0 : catch (AipsError &x)
351 : {
352 0 : logIO() << LogOrigin("PointingCal","SDS::solve")
353 : << x.getMesg()
354 0 : << LogIO::POST;
355 0 : return -1;
356 : }
357 0 : Chisq0 = Chisq;
358 : }
359 : else
360 : {
361 0 : oldOffsets = epj.solveRPar();
362 : }
363 : }
364 : // logIO() << LogOrigin("PointingCal","SDS::solve")
365 : // << "Solution interval = " << SlotNo << " Final Chisq = " << Chisq
366 : // << LogIO::NORMAL3 << LogIO::POST;
367 :
368 : //
369 : // Finalize the solutions in VisJones internal cache.
370 : //
371 0 : Chisq = fabs(Chisq0-Chisq);
372 :
373 : logIO() << LogIO::NORMAL3 << "No. of iterations used: "
374 : << iter
375 : << " Time total, per iteration: "
376 0 : << timer.all() << ", " << timer.all()/(iter+1) << " sec" << LogIO::POST;
377 :
378 :
379 0 : return (Chisq < tolerance()? Chisq:-Chisq);
380 : }
381 :
382 0 : Double SteepestDescentSolver::solve2(VisEquation& ve, VisIter& vi, EPJones& epj,
383 : Int nAnt, Int SlotNo)
384 : {
385 0 : Cube<Float> bestSolution;
386 0 : Vector<Complex> ResidualVj, dAzVj, dElVj;
387 0 : Vector<Bool> antFlagged;
388 0 : Double Chisq = std::numeric_limits<float>::max( );
389 : Double dChisq, bestChisq, sumWt, /*Time,*/ AzHDiag, ElHDiag,localGain, coVar1, coVar2,
390 : stepSize0,stepSize1;
391 0 : Int iPol=0, nPol, iter, axis0=0,axis1=0, iCorr;
392 0 : Bool Converged=false;
393 0 : Timer timer;
394 :
395 0 : nPol=polMap_p.nelements();
396 :
397 0 : ResidualVj.resize(nAnt);
398 0 : dAzVj.resize(nAnt);
399 0 : dElVj.resize(nAnt);
400 :
401 0 : epj.guessPar();
402 0 : for(Int ia=0;ia<nAnt;ia++)
403 : {
404 :
405 : // epj.solveRPar()(0,0,ia) = -0.001;
406 : // epj.solveRPar()(2,0,ia) = 0.001;
407 : // epj.solveRPar()(1,0,ia) = 0.001;
408 : // epj.solveRPar()(3,0,ia) = -0.001;
409 : //
410 : // Solved values for the squint
411 : //
412 : // epj.solveRPar()(0,0,ia) = -0.000139;
413 : // epj.solveRPar()(2,0,ia) = 5.25e-6;
414 : // epj.solveRPar()(1,0,ia) = 4.96e-5;
415 : // epj.solveRPar()(3,0,ia) = -0.000139;
416 : }
417 0 : for(Int ip=0;ip<nPol;ip++)
418 : {
419 0 : Converged=false;
420 0 : localGain = gain();
421 :
422 0 : iPol=polMap_p[ip];
423 0 : if (iPol >= 0)
424 : {
425 0 : iCorr=iPol;
426 0 : iCorr=iPol=ip;
427 : //@#$%@#$%@#%@#%@#%@#%@#%
428 : // iPol = ip;
429 : //@#$%@#$%@#%@#%@#%@#%@#%
430 :
431 0 : epj.diffResiduals(vi,ve,residual_p, gradient0_p, gradient1_p, flags);
432 0 : Chisq = getGOF(residual_p,iCorr,sumWt,"Init");
433 0 : if (sumWt == 0) return -Chisq;
434 :
435 0 : bestChisq = Chisq;
436 0 : bestSolution.resize(epj.solveRPar().shape());
437 0 : bestSolution = epj.solveRPar();
438 : //epj.guessPar();
439 : //Time = getCurrentTimeStamp(residual_p);
440 0 : timer.mark();
441 0 : antFlagged.resize(nAnt);
442 0 : logIO() << LogOrigin("PointingCal","SDS::solve2")
443 : << "Solution interval = " << SlotNo
444 : << " Initial Chisq = " << Chisq
445 : << " Polarization = " << iPol
446 : << " SumOfWt = " << sumWt
447 0 : << LogIO::NORMAL << LogIO::POST;
448 0 : if (iPol==0) {axis0=0;axis1=2;}
449 0 : else {axis0=1; axis1=3;}
450 0 : if (iPol==0) {axis0=1;axis1=3;}
451 0 : else {axis0=0; axis1=2;}
452 0 : iter=-1;
453 0 : for (iter=0;iter<numberIterations();iter++)
454 : {
455 : //
456 : // Update the current position
457 : //
458 0 : for(Int ant=0;ant<nAnt;ant++)
459 : {
460 : // cerr << "Init: " << ant << " "
461 : // << epj.solveRPar()(axis0,0,ant) << " "
462 : // << epj.solveRPar()(axis1,0,ant) << endl;
463 0 : ResidualVj = getVj(residual_p, nAnt,ant,iCorr,sumWt,0,0);
464 0 : dAzVj = getVj(gradient0_p,nAnt,ant,iCorr,sumWt,0,1);
465 0 : dElVj = getVj(gradient1_p,nAnt,ant,iCorr,sumWt,0,1);
466 :
467 0 : antFlagged[ant]=false;
468 0 : if (sumWt > 0)
469 : {
470 0 : stepSize0 = stepSize1 = coVar1 = coVar2 = 0;
471 0 : if (sumWt > 0)
472 : {
473 0 : coVar1 = real(innerProduct(dAzVj,dAzVj))/(sumWt*sumWt);
474 0 : coVar2 = real(innerProduct(dElVj,dElVj))/(sumWt*sumWt);
475 : }
476 0 : AzHDiag = ElHDiag = 0;
477 :
478 0 : if (coVar1 > 0)
479 : {
480 0 : AzHDiag = 1.0/coVar1;
481 0 : stepSize0 = real(innerProduct(ResidualVj,dAzVj))/sumWt;
482 : // cerr << "Ant:" << ant << " " << stepSize0 << " ";
483 0 : stepSize0 *= -2*AzHDiag*localGain;
484 : // cerr << stepSize0 << " ";
485 0 : epj.solveRPar()(axis0,0,ant) = epj.solveRPar()(axis0,0,ant) + stepSize0;
486 : // epj.solveRPar()(axis0,0,ant) = 0.001;
487 : }
488 0 : if (coVar2 > 0)
489 : {
490 0 : ElHDiag = 1.0/coVar2;
491 0 : stepSize1 = real(innerProduct(ResidualVj,dElVj))/sumWt;
492 : // cerr << stepSize1 << " ";
493 0 : stepSize1 *= -2*ElHDiag*localGain;
494 0 : epj.solveRPar()(axis1,0,ant) = epj.solveRPar()(axis1,0,ant) + stepSize1;
495 : // cerr << stepSize1 << " "
496 : // << epj.solveRPar()(axis0,0,ant) << " "
497 : // << epj.solveRPar()(axis1,0,ant) << endl;
498 : // epj.solveRPar()(axis1,0,ant) = -0.001;
499 : }
500 : /*
501 : cout << "SDS Sol Loop = "
502 : << iPol << " " << iter << " " << ant << " " << axis0 << " " << axis1 << " "
503 : << epj.solveRPar()(axis0,0,ant) << " "
504 : << epj.solveRPar()(axis1,0,ant) << " "
505 : << stepSize0 << " " << stepSize1 << " "
506 : << endl;
507 : */
508 : }
509 : else
510 : {
511 0 : antFlagged[ant]=true;
512 0 : epj.solveRPar()(axis0,0,ant)=0.0;//0.001;
513 0 : epj.solveRPar()(axis1,0,ant)=0.0;//-0.001;
514 : }
515 : //epj.solveRPar()(axis0,0,ant)=0.001;
516 : //epj.solveRPar()(axis1,0,ant)=-0.001;
517 :
518 : }
519 : //
520 : // Compute the residuals and the gradients with the
521 : // updated solutions and check for convergence.
522 : //
523 0 : epj.diffResiduals(vi,ve,residual_p, gradient0_p, gradient1_p, flags);
524 0 : Chisq = getGOF(residual_p,iCorr,sumWt,"Final");
525 : /*
526 : for(Int kk=0;kk<nAnt;kk++)
527 : cout << "SDS Sol Iter = "
528 : << iPol << " " << iter << " " << axis0 << " " << axis1 << " " << kk << " "
529 : << epj.solveRPar()(axis0,0,kk) << " "
530 : << epj.solveRPar()(axis1,0,kk)
531 : << Chisq << " " << bestChisq
532 : << endl;
533 : */
534 0 : dChisq = (bestChisq-Chisq);
535 : /*
536 : cout << iter << " " << iPol << " " << Chisq << " " << bestChisq
537 : << " " << dChisq << " " << " " << bestChisq << " " << localGain << endl;
538 : */
539 :
540 0 : if ((fabs(dChisq) < tolerance()))
541 : {
542 0 : Converged=true;
543 0 : epj.setRPar(bestSolution);
544 0 : Chisq = bestChisq;
545 0 : break;
546 : }
547 0 : if (Chisq >= bestChisq)
548 : {
549 : // Backtracking....
550 0 : localGain /= 2.0;
551 0 : epj.setRPar(bestSolution); //oldOffsets;
552 :
553 : // IPosition shp(bestSolution.shape());
554 : // for(Int ii=0;ii<shp(2);ii++)
555 : // for(Int jj=0;jj<shp(0);jj++)
556 : // cout << "BT: " << bestSolution(jj,0,ii) << " "
557 : // << epj.solveRPar()(jj,0,ii) << endl;
558 :
559 0 : epj.diffResiduals(vi,ve,residual_p, gradient0_p, gradient1_p, flags);
560 :
561 0 : Chisq = getGOF(residual_p,iCorr,sumWt,"BACK");
562 :
563 : // cout << "BT Chisq: " << Chisq << " "
564 : // << sum(bestSolution-epj.solveRPar())
565 : // << " " << residual.nRow() << endl;
566 : // if (Chisq != bestChisq)
567 : // {
568 : // IPosition ndx(3,0);
569 :
570 : // for(Int i=0;i<residual.nRow();i++)
571 : // {
572 : // ndx(2)=i;
573 : // // if (residual.flag()(0,i) == 0)
574 : // for(Int xx=0;xx<2;xx++)
575 : // {
576 : // ndx(0)=xx;
577 : // if (
578 : // (!residual.flagCube()(ndx)) &&
579 : // (!residual.flag()(xx,ndx(2))) &&
580 : // (!residual.flagRow()(ndx(2))) &&
581 : // (residual.antenna1()(ndx(2)) != residual.antenna2()(ndx(2)))
582 : // )
583 : // cout << xx<<"*";
584 : // }
585 :
586 : // cout << "SDS: Residual: " << i
587 : // << " " << residual.modelVisCube()(0,0,i)
588 : // << " " << residual.modelVisCube()(1,0,i)
589 : // << " " << residual.visCube()(0,0,i)
590 : // << " " << residual.visCube()(1,0,i)
591 : // << " " << residual.flagRow()(i)
592 : // << " " << residual.flag()(0,i)
593 : // << " " << residual.flag()(1,i)
594 : // << " " << residual.flagCube()(0,0,i)
595 : // << " " << residual.flagCube()(1,0,i)
596 : // << " " << residual.antenna1()(i)
597 : // << " " << residual.antenna2()(i)
598 : // << endl;
599 : // }
600 : // exit(0);
601 : // }
602 : }
603 : else
604 : {
605 0 : bestSolution = epj.solveRPar();
606 0 : bestChisq = Chisq;
607 : }
608 0 : if (localGain < 1e-20) {Converged = true; break;};
609 : }
610 : // epj.solveRPar() = bestSolution;
611 : // {
612 : // IPosition shp(bestSolution.shape());
613 : // cout << "-------------------------------------------------------------------------" << endl;
614 : // for(Int ii=0;ii<shp(2);ii++)
615 : // {
616 : // cout << "Sol: ";
617 : // for(Int jj=0;jj<shp(0);jj++)
618 : // cout << ii << " " << jj << " " << bestSolution(jj,0,ii);
619 : // cout << " " << endl;
620 : // }
621 : // }
622 0 : epj.diffResiduals(vi,ve,residual_p,gradient0_p,gradient1_p,flags);
623 0 : Chisq = getGOF(residual_p,iCorr,sumWt);
624 0 : logIO() << LogOrigin("PointingCal","SDS::solve2")
625 : << "Solution interval = " << SlotNo
626 : << " Final Chisq = " << Chisq
627 : << " Polarization = " << iPol
628 : << " NIter = " << iter
629 0 : << LogIO::NORMAL << LogIO::POST;
630 : logIO() << LogIO::NORMAL3 << "No. of iterations used: "
631 : << iter
632 : << " Time total, per iteration: "
633 0 : << timer.all() << ", " << timer.all()/(iter+1) << " sec" << LogIO::POST;
634 :
635 : } // if (iPol >= 0)
636 : } // for(Int ip=0; ip<nPol; ip++)
637 : // return (Chisq < tolerance()? Chisq:-Chisq);
638 0 : return (Converged ? Chisq: -Chisq);
639 : }
640 :
641 : };
|
