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