Line data Source code
1 : //# IncCEMemModel.cc: this implements IncCEMemModel
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,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 : //# $Id$
27 :
28 : #include <synthesis/MeasurementEquations/IncCEMemModel.h>
29 : #include <synthesis/MeasurementEquations/CEMemProgress.h>
30 : #include <synthesis/MeasurementEquations/LatConvEquation.h>
31 : #include <casacore/casa/Arrays/IPosition.h>
32 : #include <casacore/lattices/LEL/LatticeExpr.h>
33 : #include <casacore/lattices/LEL/LatticeExprNode.h>
34 : #include <casacore/lattices/Lattices/TiledLineStepper.h>
35 : #include <casacore/lattices/Lattices/TempLattice.h>
36 : #include <casacore/lattices/Lattices/LatticeIterator.h>
37 : #include <casacore/casa/Arrays/Slice.h>
38 : #include <casacore/casa/Arrays/Matrix.h>
39 : #include <casacore/casa/Arrays/Vector.h>
40 : #include <casacore/casa/Arrays/Array.h>
41 : #include <casacore/casa/Arrays/ArrayError.h>
42 : #include <casacore/casa/Arrays/ArrayMath.h>
43 : #include <casacore/casa/Arrays/VectorIter.h>
44 : #include <casacore/casa/Logging/LogOrigin.h>
45 : #include <casacore/casa/BasicMath/Math.h>
46 : #include <casacore/casa/Exceptions/Error.h>
47 : #include <casacore/casa/Utilities/Assert.h>
48 : #include <casacore/casa/System/PGPlotter.h>
49 : #include <casacore/coordinates/Coordinates.h>
50 : #include <iostream>
51 : #include <sstream>
52 :
53 :
54 : using namespace casacore;
55 : namespace casa { //# NAMESPACE CASA - BEGIN
56 :
57 : //----------------------------------------------------------------------
58 0 : IncCEMemModel::IncCEMemModel(IncEntropy &ent,
59 : Lattice<Float> & model,
60 : Lattice<Float> & deltaModel,
61 : uInt nIntegrations,
62 : Float sigma,
63 : Float targetFlux,
64 : Bool useFluxConstraint,
65 : Bool initializeModel,
66 0 : Bool doImageplane) :
67 : itsEntropy_ptr(&ent),
68 : itsResidualEquation_ptr(0),
69 : itsModel_ptr(&model),
70 : itsDeltaModel_ptr(&deltaModel),
71 : itsPrior_ptr(0),
72 : itsMask_ptr(0),
73 : itsStep_ptr(0),
74 : itsResidual_ptr(0),
75 : itsInitializeModel(initializeModel),
76 : itsNumberIterations(nIntegrations),
77 : itsSigma(sigma),
78 : itsTargetFlux(targetFlux),
79 : itsQ(10.0),
80 : itsUseFluxConstraint(useFluxConstraint),
81 : itsDoImagePlane(doImageplane),
82 : itsThreshold0(0.0),
83 : itsThresholdSpeedup(0.0),
84 : itsAlpha(0.0),
85 : itsBeta(0.0),
86 : itsFlux(0.0),
87 : itsModelFlux(0.0),
88 : itsDeltaFlux(0.0),
89 : itsFirstIteration(0),
90 : itsChoose(false),
91 0 : itsLog(LogOrigin("IncCEMemModel",
92 : "IncCEMemModel(const Lattice<Float> & model)")),
93 0 : itsProgressPtr(0)
94 : {
95 0 : initStuff();
96 0 : };
97 :
98 : //----------------------------------------------------------------------
99 0 : IncCEMemModel::IncCEMemModel(IncEntropy &ent,
100 : Lattice<Float> & model,
101 : Lattice<Float> & deltaModel,
102 : Lattice<Float> & prior,
103 : uInt nIntegrations,
104 : Float sigma,
105 : Float targetFlux,
106 : Bool useFluxConstraint,
107 : Bool initializeModel,
108 0 : Bool doImageplane) :
109 : itsEntropy_ptr(&ent),
110 : itsResidualEquation_ptr(0),
111 : itsModel_ptr(&model),
112 : itsDeltaModel_ptr(&deltaModel),
113 : itsPrior_ptr(&prior),
114 : itsMask_ptr(0),
115 : itsStep_ptr(0),
116 : itsResidual_ptr(0),
117 : itsInitializeModel(initializeModel),
118 : itsNumberIterations(nIntegrations),
119 : itsSigma(sigma),
120 : itsTargetFlux(targetFlux),
121 : itsQ(10.0),
122 : itsUseFluxConstraint(useFluxConstraint),
123 : itsDoImagePlane(doImageplane),
124 : itsThreshold0(0.0),
125 : itsThresholdSpeedup(0.0),
126 : itsAlpha(0.0),
127 : itsBeta(0.0),
128 : itsFlux(0.0),
129 : itsModelFlux(0.0),
130 : itsDeltaFlux(0.0),
131 : itsFirstIteration(0),
132 : itsChoose(false),
133 0 : itsLog(LogOrigin("IncCEMemModel",
134 : "IncCEMemModel(const Lattice<Float> & model)")),
135 0 : itsProgressPtr(0)
136 : {
137 0 : initStuff();
138 0 : };
139 :
140 : //----------------------------------------------------------------------
141 0 : IncCEMemModel::IncCEMemModel(IncEntropy &ent,
142 : Lattice<Float> & model,
143 : Lattice<Float> & deltaModel,
144 : Lattice<Float> & prior,
145 : Lattice<Float> & mask,
146 : uInt nIntegrations,
147 : Float sigma,
148 : Float targetFlux,
149 : Bool useFluxConstraint,
150 : Bool initializeModel,
151 0 : Bool doImageplane)
152 : :
153 : itsEntropy_ptr(&ent),
154 : itsResidualEquation_ptr(0),
155 : itsModel_ptr(&model),
156 : itsDeltaModel_ptr(&deltaModel),
157 : itsPrior_ptr(&prior),
158 : itsMask_ptr(&mask),
159 : itsStep_ptr(0),
160 : itsResidual_ptr(0),
161 : itsInitializeModel(initializeModel),
162 : itsNumberIterations(nIntegrations),
163 : itsSigma(sigma),
164 : itsTargetFlux(targetFlux),
165 : itsQ(10.0),
166 : itsUseFluxConstraint(useFluxConstraint),
167 : itsDoImagePlane(doImageplane),
168 : itsThreshold0(0.0),
169 : itsThresholdSpeedup(0.0),
170 : itsAlpha(0.0),
171 : itsBeta(0.0),
172 : itsFlux(0.0),
173 : itsModelFlux(0.0),
174 : itsDeltaFlux(0.0),
175 : itsFirstIteration(0),
176 : itsChoose(false),
177 0 : itsLog(LogOrigin("IncCEMemModel",
178 : "IncCEMemModel(const Lattice<Float> & model)")),
179 0 : itsProgressPtr(0)
180 :
181 : {
182 0 : initStuff();
183 0 : };
184 :
185 : //----------------------------------------------------------------------
186 0 : IncCEMemModel::IncCEMemModel(IncEntropy &ent,
187 : Lattice<Float> & model,
188 : Lattice<Float> & deltaModel,
189 : uInt nIntegrations,
190 : Lattice<Float> & mask,
191 : Float sigma,
192 : Float targetFlux,
193 : Bool useFluxConstraint,
194 : Bool initializeModel,
195 0 : Bool doImageplane)
196 : :
197 : itsEntropy_ptr(&ent),
198 : itsResidualEquation_ptr(0),
199 : itsModel_ptr(&model),
200 : itsDeltaModel_ptr(&deltaModel),
201 : itsPrior_ptr(0),
202 : itsMask_ptr(&mask),
203 : itsStep_ptr(0),
204 : itsResidual_ptr(0),
205 : itsInitializeModel(initializeModel),
206 : itsNumberIterations(nIntegrations),
207 : itsSigma(sigma),
208 : itsTargetFlux(targetFlux),
209 : itsQ(10.0),
210 : itsUseFluxConstraint(useFluxConstraint),
211 : itsDoImagePlane(doImageplane),
212 : itsThreshold0(0.0),
213 : itsThresholdSpeedup(0.0),
214 : itsAlpha(0.0),
215 : itsBeta(0.0),
216 : itsFlux(0.0),
217 : itsModelFlux(0.0),
218 : itsDeltaFlux(0.0),
219 : itsFirstIteration(0),
220 : itsChoose(false),
221 0 : itsLog(LogOrigin("IncCEMemModel",
222 : "IncCEMemModel(const Lattice<Float> & model)")),
223 0 : itsProgressPtr(0)
224 :
225 : {
226 0 : initStuff();
227 0 : };
228 :
229 :
230 : //----------------------------------------------------------------------
231 0 : IncCEMemModel::~IncCEMemModel()
232 : {
233 0 : delete itsStep_ptr; itsStep_ptr=0;
234 0 : delete itsResidual_ptr; itsResidual_ptr=0;
235 0 : };
236 :
237 :
238 :
239 : //----------------------------------------------------------------------
240 0 : void IncCEMemModel::setPrior(Lattice<Float> & prior)
241 : {
242 0 : checkImages(itsModel_ptr, &prior);
243 0 : itsPrior_ptr = &prior;
244 0 : initStuff();
245 0 : }
246 :
247 : //----------------------------------------------------------------------
248 0 : Float IncCEMemModel::getThreshold()
249 : {
250 0 : if (itsThresholdSpeedup<= 0.0) {
251 0 : return itsThreshold0;
252 : } else {
253 0 : Float factor = pow(2.0, (Float)(itsIteration - itsFirstIteration)
254 0 : / ((Float)(itsThresholdSpeedup)) );
255 0 : return (itsThreshold0 * factor);
256 : }
257 : };
258 :
259 :
260 : //----------------------------------------------------------------------
261 0 : Bool IncCEMemModel::initStuff()
262 : {
263 0 : checkImage(itsModel_ptr);
264 0 : checkImage(itsDeltaModel_ptr);
265 0 : checkImages(itsModel_ptr, itsDeltaModel_ptr);
266 :
267 : // initialize various parameters
268 :
269 0 : if (itsMask_ptr) {
270 0 : LatticeExprNode LEN = sum(*itsMask_ptr);
271 0 : itsNumberPixels = LEN.getFloat();
272 : } else {
273 0 : itsNumberPixels = itsModel_ptr->nelements();
274 : }
275 :
276 0 : formFlux();
277 :
278 0 : String entType;
279 0 : itsEntropy_ptr->entropyName(entType);
280 :
281 0 : Float clipVal = 1.0e-7;
282 0 : if (itsPrior_ptr != 0) {
283 0 : itsDefaultLevel = 0.0;
284 0 : if (entType == "ENTROPY") {
285 0 : Lattice<Float> &prior = *itsPrior_ptr;
286 0 : prior.copyData( (LatticeExpr<Float>) (iif(prior < clipVal, clipVal, prior)) );
287 : }
288 : } else {
289 0 : itsDefaultLevel = abs(itsTargetFlux) / itsNumberPixels;
290 : }
291 :
292 0 : Lattice<Float> &model = *itsModel_ptr;
293 0 : Lattice<Float> &delta = *itsDeltaModel_ptr;
294 0 : if (itsInitializeModel) {
295 0 : if (itsFirstIteration == 0) {
296 0 : delta.set(itsDefaultLevel);
297 0 : if (entType == "ENTROPY") {
298 0 : delta.copyData( (LatticeExpr<Float>) (iif( (model+delta)<clipVal, (clipVal-model), delta)) );
299 : }
300 : } else {
301 0 : delta.set(0.0);
302 : }
303 : } else {
304 0 : if (entType == "ENTROPY") {
305 : // use the model passed in, but clip it to avoid negativity
306 0 : delta.copyData( (LatticeExpr<Float>) (iif((model+delta)<clipVal, (clipVal-model), delta)) );
307 : }
308 : }
309 :
310 0 : applyMask( model ); // zero out masked pixels in model
311 :
312 0 : itsTolerance = 0.05;
313 0 : itsGain = 0.3;
314 0 : itsMaxNormGrad = 100.0;
315 :
316 0 : itsDoInit = true;
317 0 : Bool isOK = true;
318 :
319 : //Create temporary images
320 :
321 0 : if(itsStep_ptr) delete itsStep_ptr; itsStep_ptr=0;
322 0 : itsStep_ptr = new TempLattice<Float> (itsModel_ptr->shape(), 2);
323 0 : if(itsResidual_ptr) delete itsResidual_ptr; itsResidual_ptr=0;
324 0 : itsResidual_ptr = new TempLattice<Float> (itsModel_ptr->shape(), 2);
325 :
326 0 : if (itsStep_ptr && itsResidual_ptr) {
327 0 : itsStep_ptr->set(0.0);
328 0 : itsResidual_ptr->set(0.0);
329 : } else {
330 0 : isOK = false;
331 : }
332 :
333 : // We have been given an Entropy object, now we have to
334 : // tell it about US
335 0 : itsEntropy_ptr->setMemModel( *this );
336 :
337 0 : return isOK;
338 : };
339 :
340 : //----------------------------------------------------------------------
341 0 : void IncCEMemModel::initializeAlphaBeta()
342 : {
343 : // Note: in practice, this routine seems to never get called;
344 : // included for consistency with SDE MEM routine
345 0 : if (! itsUseFluxConstraint) {
346 0 : itsAlpha = max(0.0, itsGDG(H,C)/ itsGDG(C,C) );
347 0 : itsBeta = 0.0;
348 : } else {
349 0 : Double det = itsGDG(C,C)*itsGDG(F,F) - itsGDG(C,F)*itsGDG(C,F);
350 0 : itsAlpha = (itsGDG(F,F)*itsGDG(H,C) - itsGDG(C,F)*itsGDG(H,F) )/det;
351 0 : itsBeta = (itsGDG(C,C)*itsGDG(H,F) - itsGDG(C,F)*itsGDG(H,C) )/det;
352 : }
353 0 : };
354 :
355 : //----------------------------------------------------------------------
356 0 : void IncCEMemModel::updateAlphaBeta()
357 : {
358 :
359 : // code stolen from SDE's mem.f
360 0 : Double a = itsGDG(C,J)/itsGDG(C,C);
361 0 : Double b = square(a) - (itsGDG(J, J) - itsGain*itsLength)/ itsGDG(C, C);
362 :
363 : Double dAMax;
364 : Double dAMin;
365 : Double dBMax;
366 : Double dBMin;
367 :
368 0 : if (b > 0.0) {
369 0 : b = sqrt(b);
370 0 : dAMax = a + b;
371 0 : dAMin = a - b;
372 : } else {
373 0 : dAMax = 0.0;
374 0 : dAMin = 0.0;
375 : }
376 :
377 : Double dChisq;
378 : Double dAlpha;
379 : Double dBeta;
380 : Double dFlux;
381 :
382 0 : if ( ! itsUseFluxConstraint ) {
383 0 : dChisq = itsChisq - itsTargetChisq + itsGDG (C, J);
384 0 : dAlpha = dChisq/itsGDG(C, C);
385 :
386 0 : dAlpha = max(dAMin, min(dAMax, dAlpha));
387 0 : itsAlpha = max(0.0, itsAlpha + dAlpha);
388 : } else {
389 0 : a = itsGDG(F,J)/itsGDG(F,F);
390 0 : b = square(a) - (itsGDG(J, J) - itsGain*itsLength)/itsGDG(F,F);
391 0 : if ( b > 0.0) {
392 0 : b = sqrt((double)b);
393 0 : dBMax = a + b;
394 0 : dBMin = a - b;
395 : } else {
396 0 : dBMax = 0.0;
397 0 : dBMin = 0.0;
398 : }
399 :
400 0 : dChisq = itsChisq - itsTargetChisq + itsGDG(C, J);
401 0 : dFlux = itsFlux - itsTargetFlux + itsGDG(F, J);
402 0 : Double det = itsGDG(C, C)*itsGDG(F, F) - square(itsGDG(F, C));
403 0 : dAlpha = (itsGDG(F, F)*dChisq - itsGDG(C, F)*dFlux)/det;
404 0 : dBeta = (itsGDG(C, C)*dFlux - itsGDG(C, F)*dChisq)/det;
405 :
406 0 : dAlpha = max(dAMin, min(dAMax, dAlpha));
407 0 : itsAlpha = max(0.0, itsAlpha + dAlpha);
408 0 : dBeta = max(dBMin, min(dBMax, dBeta));
409 0 : itsBeta = itsBeta + dBeta;
410 : }
411 :
412 0 : };
413 :
414 :
415 :
416 : //----------------------------------------------------------------------
417 0 : void IncCEMemModel::changeAlphaBeta()
418 : {
419 0 : formGDG();
420 0 : itsLength = itsGDG(H,H) + square(itsAlpha) * itsGDG(C,C)
421 0 : + square(itsBeta) * itsGDG(F,F);
422 0 : if (itsAlpha == 0.0 && itsBeta == 0.0) {
423 0 : itsLength = itsGDG(F,F);
424 : }
425 0 : itsNormGrad = itsGDG(J,J) / itsLength;
426 0 : if (itsAlpha == 0.0)
427 0 : itsNormGrad = 0.0;
428 0 : if (itsNormGrad < itsGain) {
429 0 : updateAlphaBeta();
430 : } else {
431 0 : initializeAlphaBeta();
432 : }
433 0 : };
434 :
435 :
436 :
437 : //----------------------------------------------------------------------
438 0 : Bool IncCEMemModel::checkImage(const Lattice<Float> * /*im*/)
439 : {
440 : // I guess we don't have anything to do
441 0 : return true;
442 : };
443 :
444 :
445 :
446 : //----------------------------------------------------------------------
447 0 : Bool IncCEMemModel::checkImages(const Lattice<Float> *one, const Lattice<Float> *other)
448 : {
449 0 : Bool isOK = true;
450 :
451 0 : for (uInt i = 0; i < one->ndim(); i++) {
452 0 : AlwaysAssert(one->shape()(i) == other->shape()(i), AipsError);
453 : }
454 0 : return isOK;
455 : };
456 :
457 :
458 : //----------------------------------------------------------------------
459 0 : Bool IncCEMemModel::ok()
460 : {
461 0 : Bool isOK = true;
462 0 : if (!itsModel_ptr) {
463 0 : isOK = false;
464 0 : } else if (!itsDeltaModel_ptr) {
465 0 : isOK = false;
466 : } else {
467 0 : isOK = (checkImage(itsModel_ptr) && checkImage(itsDeltaModel_ptr));
468 0 : checkImages(itsModel_ptr, itsDeltaModel_ptr);
469 0 : if (itsPrior_ptr) {
470 0 : checkImages(itsModel_ptr, itsPrior_ptr);
471 : }
472 0 : if (itsMask_ptr) {
473 0 : checkImages(itsModel_ptr, itsMask_ptr);
474 : }
475 0 : if (itsStep_ptr) {
476 0 : checkImages(itsModel_ptr, itsStep_ptr);
477 : }
478 0 : if (itsResidual_ptr) {
479 0 : checkImages(itsModel_ptr, itsResidual_ptr);
480 : }
481 0 : if (! itsEntropy_ptr) {
482 0 : isOK = false;
483 : }
484 :
485 : // Also need to check state variables in the future!
486 : }
487 0 : return isOK;
488 : };
489 :
490 0 : void IncCEMemModel::setMask(Lattice<Float> & mask)
491 : {
492 0 : checkImages(itsModel_ptr, &mask);
493 0 : itsMask_ptr = &mask;
494 0 : initStuff();
495 0 : }
496 :
497 :
498 0 : void IncCEMemModel::state()
499 : {
500 0 : if (itsPrior_ptr == 0) {
501 0 : itsLog << "Using blank prior image" << LogIO::POST;
502 : } else {
503 0 : itsLog << "Using prior image " << LogIO::POST;
504 : }
505 0 : if (itsMask_ptr != 0) {
506 0 : itsLog << "Using mask to restrict emission" << LogIO::POST;
507 : }
508 0 : }
509 :
510 : //----------------------------------------------------------------------
511 0 : Bool IncCEMemModel::solve(ResidualEquation<Lattice<Float> > & eqn)
512 : {
513 :
514 0 : itsResidualEquation_ptr = &eqn;
515 0 : Bool converged = false;
516 0 : Bool endNow = false;
517 0 : state();
518 :
519 0 : itsEntropy_ptr->infoBanner();
520 :
521 0 : for (itsIteration = itsFirstIteration;
522 0 : ( (itsIteration < itsNumberIterations) && !endNow);
523 0 : itsIteration++) {
524 0 : if (itsDoImagePlane) {
525 0 : itsTargetChisq = square(itsSigma) * itsNumberPixels;
526 : } else {
527 0 : itsTargetChisq = square(itsSigma) * itsNumberPixels / itsQ;
528 : }
529 :
530 0 : oneIteration();
531 :
532 0 : itsFit = sqrt((double)( itsChisq / itsTargetChisq) );
533 0 : itsAFit = max(itsFit, 1.0f) *
534 0 : sqrt((double)( itsTargetChisq / itsNumberPixels) );
535 0 : itsEntropy_ptr->infoPerIteration(itsIteration);
536 :
537 0 : converged = itsEntropy_ptr->testConvergence();
538 0 : if (! converged && getThreshold() > 0.0) {
539 0 : converged = testConvergenceThreshold();
540 : }
541 :
542 0 : if (itsNormGrad > itsMaxNormGrad) {
543 0 : endNow = true;
544 0 : itsLog << " Excessive gradient: stopping now" << LogIO::EXCEPTION;
545 : }
546 :
547 0 : if (converged) {
548 0 : endNow = true;;
549 0 : itsLog << "Converged at iteration " << itsIteration+1 << LogIO::POST;
550 : }
551 : }
552 0 : if (!converged) {
553 0 : itsLog << "MEM failed to converge after " << itsNumberIterations+1
554 0 : << " iterations" << LogIO::POST;
555 : }
556 :
557 0 : formFlux();
558 0 : return converged;
559 : };
560 :
561 0 : Bool IncCEMemModel::testConvergenceThreshold()
562 : {
563 0 : Bool less = false;
564 0 : if (getThreshold() > 0.0) {
565 0 : if (itsCurrentPeakResidual < getThreshold() ) {
566 0 : less = true;
567 : }
568 : }
569 0 : return less;
570 : };
571 :
572 0 : Bool IncCEMemModel::applyMask( Lattice<Float> & lat ) {
573 0 : if (itsMask_ptr) {
574 0 : LatticeExpr<Float> exp = ( (*itsMask_ptr) * (lat) );
575 0 : lat.copyData( exp );
576 0 : return true;
577 : } else {
578 0 : return false;
579 : }
580 : };
581 :
582 :
583 0 : void IncCEMemModel::oneIteration()
584 : {
585 0 : ok();
586 :
587 0 : if (itsDoInit) {
588 0 : itsDoInit = false;
589 : // passing *this reverts to the LinearModel from which we are derived
590 0 : if (itsMask_ptr) {
591 0 : itsResidualEquation_ptr->residual( *itsResidual_ptr,
592 0 : itsChisq,
593 0 : *itsMask_ptr,
594 0 : *this);
595 : } else {
596 0 : itsResidualEquation_ptr->residual( *itsResidual_ptr,
597 0 : itsChisq,
598 0 : *this);
599 : }
600 0 : if (!itsDoImagePlane) itsChisq /= itsQ;
601 :
602 0 : if (itsAlpha == 0.0 )
603 0 : changeAlphaBeta();
604 : }
605 :
606 0 : calculateStep();
607 0 : relaxMin();
608 :
609 : // limit the step to less than the gain
610 0 : Float scale = 1.0;
611 0 : Float scalem = 1.0;
612 0 : if (itsNormGrad > 0.0)
613 0 : scalem = itsGain/itsNormGrad;
614 0 : scale = min(1.0f, scalem);
615 :
616 0 : takeStep(1.0, scale);
617 :
618 : // passing *this reverts to the LinearModel from which we are derived
619 0 : if (itsMask_ptr) {
620 0 : itsResidualEquation_ptr->residual( *itsResidual_ptr,
621 0 : itsChisq,
622 0 : *itsMask_ptr,
623 0 : *this);
624 : } else {
625 0 : itsResidualEquation_ptr->residual( *itsResidual_ptr,
626 0 : itsChisq,
627 0 : *this);
628 : }
629 0 : if (!itsDoImagePlane) itsChisq /= itsQ;
630 :
631 0 : if (itsThreshold0 > 0.0) {
632 0 : LatticeExprNode LEN = max(*itsResidual_ptr);
633 0 : Float rmax = LEN.getFloat();
634 0 : LEN = min(*itsResidual_ptr);
635 0 : Float rmin = LEN.getFloat();
636 0 : itsCurrentPeakResidual = max (abs(rmax), abs(rmin));
637 : }
638 :
639 : // form Gradient dot Step
640 0 : formGDS();
641 :
642 :
643 : // determine optimal step
644 0 : Double eps = 1.0;
645 0 : if (itsGradDotStep0 != itsGradDotStep1)
646 0 : eps = itsGradDotStep0/(itsGradDotStep0-itsGradDotStep1);
647 0 : if (scale != 0.0f) eps = min(eps, (Double)(scalem/scale) );
648 0 : if (eps <= 0.0) {
649 0 : eps = 1.0;
650 0 : itsLog << LogIO::WARN << "Step is zero" << LogIO::POST;
651 : }
652 :
653 : // Step to optimum point
654 :
655 0 : if (abs(eps-1.0) > itsGain) {
656 0 : takeStep(1.0, scale*(eps-1.0));
657 :
658 : // passing *this reverts to the LinearModel from which we are derived
659 0 : if (itsMask_ptr) {
660 0 : itsResidualEquation_ptr->residual( *itsResidual_ptr,
661 0 : itsChisq,
662 0 : *itsMask_ptr,
663 0 : *this);
664 : } else {
665 0 : itsResidualEquation_ptr->residual( *itsResidual_ptr,
666 0 : itsChisq,
667 0 : *this);
668 : }
669 :
670 0 : if (!itsDoImagePlane) itsChisq /= itsQ;
671 :
672 0 : if (itsThreshold0 > 0.0) {
673 0 : LatticeExprNode LEN = max(*itsResidual_ptr);
674 0 : Float rmax = LEN.getFloat();
675 0 : LEN = min(*itsResidual_ptr);
676 0 : Float rmin = LEN.getFloat();
677 0 : itsCurrentPeakResidual = max (abs(rmax), abs(rmin));
678 : }
679 : }
680 :
681 0 : formEntropy();
682 0 : formFlux();
683 :
684 : // readjust beam volume
685 0 : itsQ = itsQ*(1.0/max(0.5, min(2.0,eps))+1.0)/2.0;
686 :
687 0 : changeAlphaBeta();
688 0 : };
689 :
690 :
691 :
692 : //----------------------------------------------------------------------
693 0 : void IncCEMemModel::calculateStep()
694 : {
695 :
696 0 : formGDGStep();
697 :
698 0 : formFlux();
699 0 : itsGradDotStep0 = itsGDG(J,J);
700 0 : itsLength = itsGDG(H,H) + square(itsAlpha) * itsGDG(C,C)
701 0 : + square(itsBeta) * itsGDG(F,F);
702 0 : if (itsLength <= 0.0) {
703 0 : itsLength = itsGDG(F,F);
704 : }
705 0 : itsNormGrad = itsGDG(J,J) / itsLength;
706 0 : };
707 :
708 :
709 : //----------------------------------------------------------------------
710 0 : void IncCEMemModel::takeStep(Float wt1, Float wt2)
711 : {
712 0 : Lattice<Float> &model = *itsModel_ptr;
713 0 : Lattice<Float> &delta = *itsDeltaModel_ptr;
714 0 : Lattice<Float> &step = *itsStep_ptr;
715 0 : LatticeExprNode node;
716 0 : String entType;
717 0 : itsEntropy_ptr->entropyName(entType);
718 0 : if (entType == "ENTROPY") {
719 :
720 0 : itsRequiredModelMin = max (itsRequiredModelMin, 1.0e-7);
721 :
722 0 : node = iif( (model + wt1*delta + wt2*step) > itsRequiredModelMin,
723 0 : (wt1*delta + wt2*step),
724 0 : itsRequiredModelMin - model);
725 0 : delta.copyData( (LatticeExpr<Float>)node);
726 0 : } else if (entType == "EMPTINESS") {
727 :
728 0 : delta.copyData( (LatticeExpr<Float>)( wt1*delta + wt2*step ) );
729 :
730 : } else {
731 0 : cout << "UNKNOWN ENTROPY TYPE" << endl;
732 0 : AlwaysAssert( ((entType == "ENTROPY") || (entType == "EMPTINESS")),
733 : AipsError);
734 : }
735 :
736 0 : applyMask( delta );
737 0 : };
738 :
739 :
740 :
741 : //----------------------------------------------------------------------
742 0 : Float IncCEMemModel::formFlux()
743 : {
744 0 : itsFlux = 0.0;
745 0 : itsModelFlux = 0.0;
746 0 : itsDeltaFlux = 0.0;
747 0 : itsModelMin = 1e+20;
748 0 : itsModelMax = -1e+20;
749 :
750 0 : Lattice<Float> &model = *itsModel_ptr;
751 0 : Lattice<Float> &delta = *itsDeltaModel_ptr;
752 0 : TiledLineStepper mtls(model.shape(), model.niceCursorShape(), 0);
753 0 : RO_LatticeIterator<Float> mod( model, mtls);
754 0 : RO_LatticeIterator<Float> del( delta, mtls);
755 0 : Float pixsum = 0.0;
756 :
757 0 : for (mod.reset(), del.reset(); !mod.atEnd(); mod++, del++) {
758 0 : for (uInt i=0;i<mod.vectorCursor().nelements();i++) {
759 0 : itsModelFlux += mod.vectorCursor()(i);
760 0 : itsDeltaFlux += del.vectorCursor()(i);
761 0 : pixsum = mod.vectorCursor()(i) + del.vectorCursor()(i);
762 0 : if (pixsum > itsModelMax)
763 0 : itsModelMax = pixsum;
764 0 : if (pixsum < itsModelMin)
765 0 : itsModelMin = pixsum;
766 : }
767 : }
768 0 : itsFlux = itsDeltaFlux + itsModelFlux;
769 :
770 0 : return itsFlux;
771 : };
772 :
773 :
774 :
775 : } //# NAMESPACE CASA - END
776 :
|
