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 2 : CEMemModel::CEMemModel(Entropy &ent,
58 : Lattice<Float> & model,
59 : uInt nIntegrations,
60 : Float sigma,
61 : Float targetFlux,
62 : Bool useFluxConstraint,
63 : Bool initializeModel,
64 2 : Bool doImageplane) :
65 2 : itsEntropy_ptr(&ent),
66 2 : itsResidualEquation_ptr(0),
67 2 : itsModel_ptr(&model),
68 2 : itsPrior_ptr(0),
69 2 : itsMask_ptr(0),
70 2 : itsInitializeModel(initializeModel),
71 2 : itsNumberIterations(nIntegrations),
72 2 : itsSigma(sigma),
73 2 : itsTargetFlux(targetFlux),
74 2 : itsUseFluxConstraint(useFluxConstraint),
75 2 : itsDoImagePlane(doImageplane),
76 2 : itsThreshold0(0.0),
77 2 : itsThresholdSpeedup(0.0),
78 2 : itsCycleFlux(0.0),
79 2 : itsFirstIteration(0),
80 2 : itsChoose(false),
81 2 : itsLog(LogOrigin("CEMemModel",
82 : "CEMemModel(const Lattice<Float> & model)")),
83 2 : itsProgressPtr(0)
84 : {
85 2 : initStuff();
86 2 : };
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 2 : CEMemModel::~CEMemModel()
160 : {
161 2 : delete itsStep_ptr;
162 2 : delete itsResidual_ptr;
163 2 : };
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 9 : Float CEMemModel::getThreshold()
177 : {
178 9 : if (itsThresholdSpeedup<= 0.0) {
179 9 : 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 4 : Bool CEMemModel::initStuff()
189 : {
190 4 : checkImage(itsModel_ptr);
191 :
192 : // initialize various parameters
193 :
194 4 : if (itsMask_ptr) {
195 2 : LatticeExprNode LEN = sum(*itsMask_ptr);
196 2 : itsNumberPixels = LEN.getFloat();
197 2 : } else {
198 2 : itsNumberPixels = itsModel_ptr->nelements();
199 : }
200 :
201 4 : formFlux();
202 4 : itsDefaultLevel = abs(itsTargetFlux) / itsNumberPixels;
203 :
204 4 : Float clipVal = 1.0e-7;
205 4 : 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 4 : Lattice<Float> &model = *itsModel_ptr;
212 4 : if (itsInitializeModel) {
213 4 : if (itsDefaultLevel > 0.0) {
214 4 : 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 4 : applyMask( model ); // zero out masked pixels in model
224 :
225 4 : itsTolerance = 0.05;
226 4 : itsQ = 10;
227 4 : itsGain = 0.3;
228 4 : itsMaxNormGrad = 100.0;
229 :
230 4 : itsDoInit = true;
231 4 : itsAlpha = 0;
232 4 : itsBeta = 0;
233 4 : Bool isOK = true;
234 :
235 : //Create temporary images
236 :
237 4 : itsStep_ptr = new TempLattice<Float> (itsModel_ptr->shape(), 2);
238 4 : itsResidual_ptr = new TempLattice<Float> (itsModel_ptr->shape(), 2);
239 4 : if (itsStep_ptr && itsResidual_ptr) {
240 4 : itsStep_ptr->set(0.0);
241 4 : 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 4 : itsEntropy_ptr->setMemModel( *this );
249 :
250 4 : return isOK;
251 : };
252 :
253 :
254 : //----------------------------------------------------------------------
255 2 : void CEMemModel::initializeAlphaBeta()
256 : {
257 : // Note: in practice, this routine seems to never get called;
258 : // included for consistency with SDE MEM routine
259 2 : if (! itsUseFluxConstraint) {
260 2 : itsAlpha = max(0.0, itsGDG(H,C)/ itsGDG(C,C) );
261 2 : 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 2 : };
268 :
269 :
270 :
271 : //----------------------------------------------------------------------
272 10 : void CEMemModel::updateAlphaBeta()
273 : {
274 :
275 : // code stolen from SDE's mem.f
276 10 : Double a = itsGDG(C,J)/itsGDG(C,C);
277 10 : 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 10 : if (b > 0.0) {
285 10 : b = sqrt(b);
286 10 : dAMax = a + b;
287 10 : 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 10 : if ( ! itsUseFluxConstraint ) {
299 10 : dChisq = itsChisq - itsTargetChisq + itsGDG (C, J);
300 10 : dAlpha = dChisq/itsGDG(C, C);
301 :
302 10 : dAlpha = max(dAMin, min(dAMax, dAlpha));
303 10 : 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 10 : };
329 :
330 :
331 :
332 : //----------------------------------------------------------------------
333 12 : void CEMemModel::changeAlphaBeta()
334 : {
335 12 : formGDG();
336 12 : itsLength = itsGDG(H,H) + square(itsAlpha) * itsGDG(C,C)
337 12 : + square(itsBeta) * itsGDG(F,F);
338 12 : if (itsAlpha == 0.0 && itsBeta == 0.0) {
339 2 : itsLength = itsGDG(F,F);
340 : }
341 12 : itsNormGrad = itsGDG(J,J) / itsLength;
342 12 : if (itsAlpha == 0.0)
343 2 : itsNormGrad = 0.0;
344 12 : if (itsNormGrad < itsGain) {
345 10 : updateAlphaBeta();
346 : } else {
347 2 : initializeAlphaBeta();
348 : }
349 12 : };
350 :
351 :
352 :
353 : //----------------------------------------------------------------------
354 14 : Bool CEMemModel::checkImage(const Lattice<Float>* /*im*/)
355 : {
356 : // I guess we don't have anything to do
357 14 : return true;
358 : };
359 :
360 :
361 :
362 : //----------------------------------------------------------------------
363 32 : Bool CEMemModel::checkImages(const Lattice<Float> *one, const Lattice<Float> *other)
364 : {
365 32 : Bool isOK = true;
366 :
367 160 : for (uInt i = 0; i < one->ndim(); i++) {
368 128 : AlwaysAssert(one->shape()(i) == other->shape()(i), AipsError);
369 : }
370 32 : return isOK;
371 : };
372 :
373 :
374 : //----------------------------------------------------------------------
375 10 : Bool CEMemModel::ok()
376 : {
377 10 : Bool isOK = true;
378 10 : if (!itsModel_ptr) {
379 0 : isOK = false;
380 : } else {
381 10 : isOK = checkImage(itsModel_ptr);
382 10 : if (itsPrior_ptr) {
383 0 : checkImages(itsModel_ptr, itsPrior_ptr);
384 : }
385 10 : if (itsMask_ptr) {
386 10 : checkImages(itsModel_ptr, itsMask_ptr);
387 : }
388 10 : if (itsStep_ptr) {
389 10 : checkImages(itsModel_ptr, itsStep_ptr);
390 : }
391 10 : if (itsResidual_ptr) {
392 10 : checkImages(itsModel_ptr, itsResidual_ptr);
393 : }
394 10 : if (! itsEntropy_ptr) {
395 0 : isOK = false;
396 : }
397 :
398 : // Also need to check state variables in the future!
399 : }
400 10 : return isOK;
401 : };
402 :
403 2 : void CEMemModel::setMask(Lattice<Float> & mask)
404 : {
405 2 : checkImages(itsModel_ptr, &mask);
406 2 : itsMask_ptr = &mask;
407 2 : initStuff();
408 2 : }
409 :
410 :
411 2 : void CEMemModel::state()
412 : {
413 2 : if (itsPrior_ptr == 0) {
414 2 : itsLog << "Using blank prior image" << LogIO::POST;
415 : } else {
416 0 : itsLog << "Using prior image " << LogIO::POST;
417 : }
418 2 : if (itsMask_ptr != 0) {
419 2 : itsLog << "Using mask to restrict emission" << LogIO::POST;
420 : }
421 2 : }
422 :
423 : //----------------------------------------------------------------------
424 2 : Bool CEMemModel::solve(ResidualEquation<Lattice<Float> > & eqn)
425 : {
426 2 : itsResidualEquation_ptr = &eqn;
427 2 : Bool converged = false;
428 2 : Bool endNow = false;
429 2 : state();
430 :
431 2 : itsEntropy_ptr->infoBanner();
432 :
433 2 : for (itsIteration = itsFirstIteration;
434 12 : ( (itsIteration < itsNumberIterations) && !endNow);
435 10 : itsIteration++) {
436 10 : if (itsDoImagePlane) {
437 10 : itsTargetChisq = square(itsSigma) * itsNumberPixels;
438 : } else {
439 0 : itsTargetChisq = square(itsSigma) * itsNumberPixels / itsQ;
440 : }
441 :
442 10 : oneIteration();
443 :
444 10 : itsFit = sqrt((double)( itsChisq / itsTargetChisq) );
445 10 : itsAFit = max(itsFit, 1.0f) *
446 10 : sqrt((double)( itsTargetChisq / itsNumberPixels) );
447 10 : itsEntropy_ptr->infoPerIteration(itsIteration);
448 :
449 10 : converged = itsEntropy_ptr->testConvergence();
450 10 : if (! converged && getThreshold() > 0.0) {
451 0 : converged = testConvergenceThreshold();
452 : }
453 :
454 10 : if (itsNormGrad > itsMaxNormGrad) {
455 0 : endNow = true;
456 0 : itsLog << " Excessive gradient: stopping now" << LogIO::EXCEPTION;
457 : }
458 :
459 10 : if (converged) {
460 1 : endNow = true;;
461 1 : itsLog << "Converged at iteration " << itsIteration+1 << LogIO::POST;
462 : }
463 :
464 : }
465 2 : if (!converged) {
466 2 : itsLog << "MEM failed to converge after " << itsNumberIterations+1
467 1 : << " iterations" << LogIO::POST;
468 : }
469 :
470 2 : formFlux();
471 2 : 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 34 : Bool CEMemModel::applyMask( Lattice<Float> & lat ) {
486 34 : if (itsMask_ptr) {
487 64 : LatticeExpr<Float> exp = ( (*itsMask_ptr) * (lat) );
488 32 : lat.copyData( exp );
489 32 : return true;
490 32 : } else {
491 2 : return false;
492 : }
493 : };
494 :
495 :
496 10 : void CEMemModel::oneIteration()
497 : {
498 10 : ok();
499 10 : if (itsDoInit) {
500 2 : itsDoInit = false;
501 : // passing *this reverts to the LinearModel from which we are derived
502 2 : itsResidualEquation_ptr->residual( *itsResidual_ptr,
503 2 : itsChisq,
504 : *this);
505 2 : applyMask( *itsResidual_ptr );
506 :
507 2 : if (!itsDoImagePlane) itsChisq /= itsQ;
508 :
509 2 : if (itsAlpha == 0.0 || itsBeta == 0.0)
510 2 : changeAlphaBeta();
511 : }
512 :
513 :
514 10 : calculateStep();
515 10 : relaxMin();
516 :
517 : // limit the step to less than the gain
518 10 : Float scale = 1.0;
519 10 : Float scalem = 1.0;
520 10 : if (itsNormGrad > 0.0)
521 10 : scalem = itsGain/itsNormGrad;
522 10 : scale = min(1.0f, scalem);
523 :
524 10 : takeStep(1.0, scale);
525 :
526 : // passing *this reverts to the LinearModel from which we are derived
527 10 : itsResidualEquation_ptr->residual( *itsResidual_ptr,
528 10 : itsChisq,
529 : *this);
530 10 : applyMask( *itsResidual_ptr );
531 10 : if (!itsDoImagePlane) itsChisq /= itsQ;
532 :
533 10 : 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 10 : formGDS();
543 :
544 :
545 : // determine optimal step
546 10 : Double eps = 1.0;
547 10 : if (itsGradDotStep0 != itsGradDotStep1)
548 10 : eps = itsGradDotStep0/(itsGradDotStep0-itsGradDotStep1);
549 10 : if (scale != 0.0f) eps = min(eps, (Double)(scalem/scale) );
550 10 : 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 10 : if (abs(eps-1.0) > itsGain) {
558 4 : takeStep(1.0, scale*(eps-1.0));
559 : // passing *this reverts to the LinearModel from which we are derived
560 4 : itsResidualEquation_ptr->residual( *itsResidual_ptr,
561 4 : itsChisq,
562 : *this);
563 4 : applyMask( *itsResidual_ptr );
564 4 : if (!itsDoImagePlane) itsChisq /= itsQ;
565 :
566 4 : 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 10 : formEntropy();
577 : // formFlux();
578 :
579 : // readjust beam volume
580 10 : itsQ = itsQ*(1.0/max(0.5, min(2.0,eps))+1.0)/2.0;
581 :
582 10 : changeAlphaBeta();
583 10 : };
584 :
585 :
586 :
587 : //----------------------------------------------------------------------
588 10 : void CEMemModel::calculateStep()
589 : {
590 :
591 10 : formGDGStep();
592 :
593 10 : formFlux();
594 10 : itsGradDotStep0 = itsGDG(J,J);
595 10 : itsLength = itsGDG(H,H) + square(itsAlpha) * itsGDG(C,C)
596 10 : + square(itsBeta) * itsGDG(F,F);
597 10 : if (itsLength <= 0.0) {
598 0 : itsLength = itsGDG(F,F);
599 : }
600 10 : itsNormGrad = itsGDG(J,J) / itsLength;
601 10 : };
602 :
603 :
604 : //----------------------------------------------------------------------
605 14 : void CEMemModel::takeStep(Float wt1, Float wt2)
606 : {
607 14 : Lattice<Float> &model = *itsModel_ptr;
608 14 : Lattice<Float> &step = *itsStep_ptr;
609 14 : LatticeExprNode node;
610 :
611 14 : applyMask( step );
612 :
613 14 : node = max ( (wt1*model + wt2*step), itsRequiredModelMin);
614 14 : model.copyData( (LatticeExpr<Float>)node);
615 14 : };
616 :
617 :
618 :
619 : //----------------------------------------------------------------------
620 16 : Float CEMemModel::formFlux()
621 : {
622 16 : itsFlux = 0.0;
623 16 : itsModelMin = 1e+20;
624 16 : itsModelMax = -1e+20;
625 :
626 16 : Lattice<Float> &model = *itsModel_ptr;
627 32 : TiledLineStepper mtls(model.shape(), model.niceCursorShape(), 0);
628 16 : RO_LatticeIterator<Float> mod( model, mtls);
629 3216 : for (mod.reset(); !mod.atEnd(); mod++) {
630 643200 : for (uInt i=0;i<mod.vectorCursor().nelements();i++) {
631 640000 : itsFlux += mod.vectorCursor()(i);
632 640000 : if (mod.vectorCursor()(i) > itsModelMax)
633 2788 : itsModelMax = mod.vectorCursor()(i);
634 640000 : if (mod.vectorCursor()(i) < itsModelMin)
635 3938 : itsModelMin = mod.vectorCursor()(i);
636 : }
637 : }
638 16 : return itsFlux;
639 16 : };
640 :
641 :
642 :
643 : } //# NAMESPACE CASA - END
644 :
|
