Line data Source code
1 : //# IncEntropy.cc: this implements IncEntropy
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2002
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/IncCEMemModel.h>
29 : #include <synthesis/MeasurementEquations/CEMemProgress.h>
30 : #include <synthesis/MeasurementEquations/IncEntropy.h>
31 : #include <casacore/lattices/LEL/LatticeExpr.h>
32 : #include <casacore/lattices/LEL/LatticeExprNode.h>
33 : #include <casacore/lattices/Lattices/LatticeStepper.h>
34 : #include <casacore/lattices/Lattices/LatticeIterator.h>
35 : #include <casacore/casa/BasicMath/Math.h>
36 : #include <iostream>
37 :
38 : using namespace casacore;
39 : namespace casa { //# NAMESPACE CASA - BEGIN
40 :
41 : //----------------------------------------------------------------------
42 0 : IncEntropy::IncEntropy()
43 : {
44 0 : cemem_ptr = 0;
45 0 : };
46 :
47 : //----------------------------------------------------------------------
48 0 : IncEntropy::~IncEntropy()
49 : {
50 : // cemem_ptr->letEntropyDie();
51 0 : };
52 :
53 : //----------------------------------------------------------------------
54 0 : IncEntropyI::IncEntropyI()
55 : {
56 0 : };
57 :
58 : //----------------------------------------------------------------------
59 0 : IncEntropyI::~IncEntropyI()
60 : {
61 0 : };
62 :
63 : //----------------------------------------------------------------------
64 :
65 0 : Float IncEntropyI::formEntropy()
66 : {
67 0 : Lattice<Float> &model = *(cemem_ptr->itsModel_ptr);
68 0 : Lattice<Float> &delta = *(cemem_ptr->itsDeltaModel_ptr);
69 0 : Float flux = 0.0;
70 0 : if (cemem_ptr->itsMask_ptr) {
71 0 : Lattice<Float> &mask = *(cemem_ptr->itsMask_ptr);
72 0 : flux = sum((mask * (model+delta))).getFloat();
73 : } else {
74 0 : flux = sum(model+delta).getFloat();
75 : }
76 :
77 0 : cemem_ptr->itsFlux = flux;
78 0 : Float defLev = cemem_ptr->itsDefaultLevel;
79 0 : LatticeExprNode myEntropyLEN;
80 0 : Float myEntropy = 0;
81 :
82 0 : if (cemem_ptr->itsMask_ptr) { // Mask
83 0 : Lattice<Float> &mask = *(cemem_ptr->itsMask_ptr);
84 0 : if (cemem_ptr->itsPrior_ptr) { // Mask and Prior
85 0 : Lattice<Float> &prior = *(cemem_ptr->itsPrior_ptr);
86 0 : myEntropyLEN = -sum(((model+delta) * log((model+delta)/prior))[mask>0.0] );
87 0 : myEntropy = myEntropyLEN.getFloat();
88 : } else { // Mask, No Prior
89 0 : myEntropyLEN = - sum( ((model+delta) * log(model+delta))[mask>0.0]) ;
90 0 : myEntropy = myEntropyLEN.getFloat() + flux * log(defLev);
91 : }
92 0 : if (flux > 0.0) {
93 0 : myEntropy = myEntropy/flux +
94 0 : log(cemem_ptr->itsNumberPixels);
95 : } else {
96 0 : myEntropy = 0.0;
97 : }
98 : } else { // No Mask; Oh for the joy of LEL!
99 0 : if (cemem_ptr->itsPrior_ptr) { // No Mask, But Prior
100 0 : Lattice<Float> &prior = *(cemem_ptr->itsPrior_ptr);
101 0 : myEntropyLEN = - sum( (model+delta) * log((model+delta)/prior)) ;
102 0 : myEntropy = myEntropyLEN.getFloat();
103 : } else { // No Mask, No Prior
104 0 : myEntropyLEN = - sum( (model+delta) * log(model+delta)) ;
105 0 : myEntropy = myEntropyLEN.getFloat() + flux * log(defLev);
106 : }
107 0 : if (flux > 0.0) {
108 0 : myEntropy = myEntropy/flux +
109 0 : log(cemem_ptr->itsNumberPixels);
110 : } else {
111 0 : myEntropy = 0.0;
112 : }
113 : }
114 :
115 0 : return myEntropy;
116 0 : };
117 :
118 :
119 : //----------------------------------------------------------------------
120 0 : void IncEntropyI::formGDG(Matrix<Double>& GDG)
121 : {
122 0 : if (cemem_ptr->itsMask_ptr) {
123 0 : Lattice<Float> &model = *(cemem_ptr->itsModel_ptr);
124 0 : Lattice<Float> &delta = *(cemem_ptr->itsDeltaModel_ptr);
125 0 : Lattice<Float> &resid = *(cemem_ptr->itsResidual_ptr);
126 0 : Lattice<Float> &mask = *(cemem_ptr->itsMask_ptr);
127 :
128 0 : GDG.resize(4,4);
129 0 : GDG.set(0.0);
130 0 : Float alpha = cemem_ptr->itsAlpha;
131 0 : Float beta = cemem_ptr->itsBeta;
132 0 : Float defLev = cemem_ptr->itsDefaultLevel;
133 :
134 : // should not use Lattice Expression Language for efficiency
135 : // using it right now to get up and running
136 0 : Float ggc = 2 * alpha * cemem_ptr->itsQ;
137 : Float rHess;
138 :
139 : RO_LatticeIterator<Float>
140 0 : mod( model, LatticeStepper(model.shape(), model.niceCursorShape()));
141 : RO_LatticeIterator<Float>
142 0 : del( delta, LatticeStepper(model.shape(), model.niceCursorShape()));
143 : RO_LatticeIterator<Float>
144 0 : res( resid, LatticeStepper(model.shape(), model.niceCursorShape()));
145 : RO_LatticeIterator<Float>
146 0 : mas( mask, LatticeStepper(model.shape(), model.niceCursorShape()));
147 : uInt i;
148 : Bool modDeleteIt;
149 : Bool delDeleteIt;
150 : Bool resDeleteIt;
151 : Bool masDeleteIt;
152 0 : Double GDGHH = 0.0;
153 0 : Double GDGHC = 0.0;
154 0 : Double GDGHF = 0.0;
155 0 : Double GDGCC = 0.0;
156 0 : Double GDGCF = 0.0;
157 0 : Double GDGFF = 0.0;
158 0 : Double gradH = 0.0;
159 0 : Double gradC = 0.0;
160 0 : if (cemem_ptr->itsPrior_ptr == 0) {
161 0 : Float logDef = log(defLev);
162 0 : for (mod.reset(), del.reset(), res.reset(), mas.reset();
163 0 : !mod.atEnd();
164 0 : mod++, del++, res++, mas++) {
165 0 : const Float *modStore = mod.cursor().getStorage(modDeleteIt);
166 0 : const Float *delStore = del.cursor().getStorage(delDeleteIt);
167 0 : const Float *resStore = res.cursor().getStorage(resDeleteIt);
168 0 : const Float *masStore = mas.cursor().getStorage(masDeleteIt);
169 0 : const Float *modIter = modStore;
170 0 : const Float *delIter = delStore;
171 0 : const Float *resIter = resStore;
172 0 : const Float *masIter = masStore;
173 0 : for (i=0;i<mod.cursor().nelements();i++,modIter++, delIter++, resIter++, masIter++) {
174 0 : if (*masIter > 0.0) {
175 0 : rHess = (*modIter + *delIter)/(1 + ggc * (*modIter + *delIter));
176 0 : gradH = -log(*modIter + *delIter) + logDef;
177 0 : gradC = -2*(*resIter);
178 0 : GDGHH = GDGHH + gradH * rHess * gradH;
179 0 : GDGHC = GDGHC + gradH * rHess * gradC;
180 0 : GDGHF = GDGHF + gradH * rHess;
181 0 : GDGCC = GDGCC + gradC * rHess * gradC;
182 0 : GDGCF = GDGCF + gradC * rHess;
183 0 : GDGFF = GDGFF + rHess;
184 : }
185 : }
186 0 : mod.cursor().freeStorage(modStore, modDeleteIt);
187 0 : del.cursor().freeStorage(delStore, delDeleteIt);
188 0 : res.cursor().freeStorage(resStore, resDeleteIt);
189 0 : mas.cursor().freeStorage(masStore, masDeleteIt);
190 : }
191 : } else {
192 0 : Lattice<Float> &prior = *(cemem_ptr->itsPrior_ptr);
193 : RO_LatticeIterator<Float>
194 0 : pri( prior, LatticeStepper(prior.shape(), model.niceCursorShape()));
195 : Bool priDeleteIt;
196 0 : for (mod.reset(), del.reset(), res.reset(), mas.reset(), pri.reset();
197 0 : !mod.atEnd();
198 0 : mod++, del++, res++, mas++, pri++) {
199 0 : const Float *modStore = mod.cursor().getStorage(modDeleteIt);
200 0 : const Float *delStore = del.cursor().getStorage(delDeleteIt);
201 0 : const Float *resStore = res.cursor().getStorage(resDeleteIt);
202 0 : const Float *masStore = mas.cursor().getStorage(masDeleteIt);
203 0 : const Float *priStore = pri.cursor().getStorage(priDeleteIt);
204 0 : const Float *modIter = modStore;
205 0 : const Float *delIter = delStore;
206 0 : const Float *resIter = resStore;
207 0 : const Float *masIter = masStore;
208 0 : const Float *priIter = priStore;
209 0 : for (i=0;i<mod.cursor().nelements();
210 0 : i++,modIter++, delIter++, resIter++, masIter++, priIter++) {
211 0 : if (*masIter > 0.0) {
212 0 : rHess = (*modIter + *delIter)/(1 + ggc* (*modIter + *delIter));
213 0 : gradH = -log( (*modIter + *delIter) / (*priIter));
214 0 : gradC = -2*(*resIter);
215 0 : GDGHH = GDGHH + gradH * rHess * gradH;
216 0 : GDGHC = GDGHC + gradH * rHess * gradC;
217 0 : GDGHF = GDGHF + gradH * rHess;
218 0 : GDGCC = GDGCC + gradC * rHess * gradC;
219 0 : GDGCF = GDGCF + gradC * rHess;
220 0 : GDGFF = GDGFF + rHess;
221 : }
222 : }
223 0 : mod.cursor().freeStorage(modStore, modDeleteIt);
224 0 : del.cursor().freeStorage(delStore, delDeleteIt);
225 0 : res.cursor().freeStorage(resStore, resDeleteIt);
226 0 : mas.cursor().freeStorage(masStore, masDeleteIt);
227 0 : pri.cursor().freeStorage(priStore, priDeleteIt);
228 : }
229 0 : }
230 0 : GDG(H,H) = GDGHH;
231 0 : GDG(H,C) = GDGHC;
232 0 : GDG(H,F) = GDGHF;
233 0 : GDG(C,C) = GDGCC;
234 0 : GDG(C,F) = GDGCF;
235 0 : GDG(F,F) = GDGFF;
236 0 : GDG(H,J) = GDGHH - alpha * GDGHC - beta * GDGHF;
237 0 : GDG(C,J) = GDGHC - alpha * GDGCC - beta * GDGCF;
238 0 : GDG(F,J) = GDGHF - alpha * GDGCF - beta * GDGFF;
239 0 : GDG(J,J) = GDGHH + square(alpha) * GDGCC
240 0 : + square(beta)*GDGFF + 2*alpha*beta*GDGCF
241 0 : - 2*alpha*GDGHC - 2*beta*GDGHF;
242 :
243 0 : } else { // no mask
244 :
245 0 : Lattice<Float> &model = *(cemem_ptr->itsModel_ptr);
246 0 : Lattice<Float> &delta = *(cemem_ptr->itsDeltaModel_ptr);
247 0 : Lattice<Float> &resid = *(cemem_ptr->itsResidual_ptr);
248 0 : GDG.resize(4,4);
249 0 : GDG.set(0.0);
250 0 : Float alpha = cemem_ptr->itsAlpha;
251 0 : Float beta = cemem_ptr->itsBeta;
252 0 : Float defLev = cemem_ptr->itsDefaultLevel;
253 :
254 : // should not use Lattice Expression Language for efficiency
255 : // using it right now to get up and running
256 0 : Float ggc = 2 * alpha * cemem_ptr->itsQ;
257 : Float rHess;
258 :
259 : RO_LatticeIterator<Float>
260 0 : mod( model, LatticeStepper(model.shape(), model.niceCursorShape()));
261 : RO_LatticeIterator<Float>
262 0 : del( delta, LatticeStepper(model.shape(), model.niceCursorShape()));
263 : RO_LatticeIterator<Float>
264 0 : res( resid, LatticeStepper(model.shape(), model.niceCursorShape()));
265 : uInt i;
266 : Bool modDeleteIt;
267 : Bool delDeleteIt;
268 : Bool resDeleteIt;
269 0 : Double GDGHH = 0.0;
270 0 : Double GDGHC = 0.0;
271 0 : Double GDGHF = 0.0;
272 0 : Double GDGCC = 0.0;
273 0 : Double GDGCF = 0.0;
274 0 : Double GDGFF = 0.0;
275 0 : Double gradH = 0.0;
276 0 : Double gradC = 0.0;
277 0 : if (cemem_ptr->itsPrior_ptr == 0) {
278 0 : Float logDef = log(defLev);
279 0 : for (mod.reset(), del.reset(), res.reset();
280 0 : !mod.atEnd();
281 0 : mod++, del++, res++) {
282 0 : const Float *modStore = mod.cursor().getStorage(modDeleteIt);
283 0 : const Float *delStore = del.cursor().getStorage(delDeleteIt);
284 0 : const Float *resStore = res.cursor().getStorage(resDeleteIt);
285 0 : const Float *modIter = modStore;
286 0 : const Float *delIter = delStore;
287 0 : const Float *resIter = resStore;
288 0 : for (i=0;i<mod.cursor().nelements();i++,modIter++, delIter++, resIter++) {
289 0 : rHess = (*modIter + *delIter)/(1 + ggc * (*modIter + *delIter));
290 0 : gradH = -log(*modIter + *delIter) + logDef;
291 0 : gradC = -2*(*resIter);
292 0 : GDGHH = GDGHH + gradH * rHess * gradH;
293 0 : GDGHC = GDGHC + gradH * rHess * gradC;
294 0 : GDGHF = GDGHF + gradH * rHess;
295 0 : GDGCC = GDGCC + gradC * rHess * gradC;
296 0 : GDGCF = GDGCF + gradC * rHess;
297 0 : GDGFF = GDGFF + rHess;
298 : }
299 0 : mod.cursor().freeStorage(modStore, modDeleteIt);
300 0 : del.cursor().freeStorage(delStore, delDeleteIt);
301 0 : res.cursor().freeStorage(resStore, resDeleteIt);
302 : }
303 : } else {
304 0 : Lattice<Float> &prior = *(cemem_ptr->itsPrior_ptr);
305 : RO_LatticeIterator<Float>
306 0 : pri( prior, LatticeStepper(prior.shape(), model.niceCursorShape()));
307 : Bool priDeleteIt;
308 0 : for (mod.reset(), del.reset(), res.reset(), pri.reset(); !mod.atEnd();
309 0 : mod++, del++, res++, pri++) {
310 0 : const Float *modStore = mod.cursor().getStorage(modDeleteIt);
311 0 : const Float *delStore = del.cursor().getStorage(delDeleteIt);
312 0 : const Float *resStore = res.cursor().getStorage(resDeleteIt);
313 0 : const Float *priStore = pri.cursor().getStorage(priDeleteIt);
314 0 : const Float *modIter = modStore;
315 0 : const Float *delIter = delStore;
316 0 : const Float *resIter = resStore;
317 0 : const Float *priIter = priStore;
318 0 : for (i=0;i<mod.cursor().nelements();i++,modIter++, delIter++, resIter++, priIter++) {
319 0 : rHess = (*modIter + *delIter)/(1 + ggc* (*modIter + *delIter));
320 0 : gradH = -log( (*modIter + *delIter) / (*priIter));
321 0 : gradC = -2*(*resIter);
322 0 : GDGHH = GDGHH + gradH * rHess * gradH;
323 0 : GDGHC = GDGHC + gradH * rHess * gradC;
324 0 : GDGHF = GDGHF + gradH * rHess;
325 0 : GDGCC = GDGCC + gradC * rHess * gradC;
326 0 : GDGCF = GDGCF + gradC * rHess;
327 0 : GDGFF = GDGFF + rHess;
328 : }
329 0 : mod.cursor().freeStorage(modStore, modDeleteIt);
330 0 : del.cursor().freeStorage(delStore, delDeleteIt);
331 0 : res.cursor().freeStorage(resStore, resDeleteIt);
332 0 : pri.cursor().freeStorage(priStore, priDeleteIt);
333 : }
334 0 : }
335 0 : GDG(H,H) = GDGHH;
336 0 : GDG(H,C) = GDGHC;
337 0 : GDG(H,F) = GDGHF;
338 0 : GDG(C,C) = GDGCC;
339 0 : GDG(C,F) = GDGCF;
340 0 : GDG(F,F) = GDGFF;
341 0 : GDG(H,J) = GDGHH - alpha * GDGHC - beta * GDGHF;
342 0 : GDG(C,J) = GDGHC - alpha * GDGCC - beta * GDGCF;
343 0 : GDG(F,J) = GDGHF - alpha * GDGCF - beta * GDGFF;
344 0 : GDG(J,J) = GDGHH + square(alpha) * GDGCC
345 0 : + square(beta)*GDGFF + 2*alpha*beta*GDGCF
346 0 : - 2*alpha*GDGHC - 2*beta*GDGHF;
347 0 : }
348 0 : };
349 :
350 :
351 :
352 : //----------------------------------------------------------------------
353 0 : void IncEntropyI::formGDGStep(Matrix<Double>& GDG)
354 : {
355 :
356 0 : if (cemem_ptr->itsMask_ptr) {
357 :
358 0 : Lattice<Float> &model = *(cemem_ptr->itsModel_ptr);
359 0 : Lattice<Float> &delta = *(cemem_ptr->itsDeltaModel_ptr);
360 0 : Lattice<Float> &resid = *(cemem_ptr->itsResidual_ptr);
361 0 : Lattice<Float> &step = *(cemem_ptr->itsStep_ptr);
362 0 : Lattice<Float> &mask = *(cemem_ptr->itsMask_ptr);
363 :
364 0 : GDG.resize(4,4);
365 0 : GDG.set(0.0);
366 0 : Float alpha = cemem_ptr->itsAlpha;
367 0 : Float beta = cemem_ptr->itsBeta;
368 0 : Float defLev = cemem_ptr->itsDefaultLevel;
369 :
370 : // should not use Lattice Expression Language for efficiency
371 : // using it right now to get up and running
372 0 : Float ggc = 2 * alpha * cemem_ptr->itsQ;
373 : Float rHess;
374 :
375 : RO_LatticeIterator<Float>
376 0 : mod( model, LatticeStepper(model.shape(), model.niceCursorShape()));
377 : RO_LatticeIterator<Float>
378 0 : del( delta, LatticeStepper(model.shape(), model.niceCursorShape()));
379 : RO_LatticeIterator<Float>
380 0 : res( resid, LatticeStepper(resid.shape(), model.niceCursorShape()));
381 : RO_LatticeIterator<Float>
382 0 : mas( mask, LatticeStepper(model.shape(), model.niceCursorShape()));
383 : LatticeIterator<Float>
384 0 : stp( step, LatticeStepper(step.shape(), model.niceCursorShape()));
385 : uInt i;
386 : Bool modDeleteIt;
387 : Bool delDeleteIt;
388 : Bool resDeleteIt;
389 : Bool stpDeleteIt;
390 : Bool masDeleteIt;
391 0 : Double GDGHH = 0.0;
392 0 : Double GDGHC = 0.0;
393 0 : Double GDGHF = 0.0;
394 0 : Double GDGCC = 0.0;
395 0 : Double GDGCF = 0.0;
396 0 : Double GDGFF = 0.0;
397 0 : Double gradH = 0.0;
398 0 : Double gradC = 0.0;
399 0 : Double gradJ = 0.0;
400 : Float stepValue;
401 0 : if ( cemem_ptr->itsPrior_ptr == 0) {
402 0 : Float logDef = log(defLev);
403 0 : for (mod.reset(), del.reset(), res.reset(), mas.reset(), stp.reset();
404 0 : !mod.atEnd();
405 0 : mod++, del++, res++, mas++, stp++) {
406 0 : const Float *modStore = mod.cursor().getStorage(modDeleteIt);
407 0 : const Float *delStore = del.cursor().getStorage(delDeleteIt);
408 0 : const Float *resStore = res.cursor().getStorage(resDeleteIt);
409 0 : const Float *masStore = mas.cursor().getStorage(masDeleteIt);
410 0 : Float *stpStore = stp.rwCursor().getStorage(stpDeleteIt);
411 0 : const Float *modIter = modStore;
412 0 : const Float *delIter = delStore;
413 0 : const Float *resIter = resStore;
414 0 : const Float *masIter = masStore;
415 0 : Float *stpIter = stpStore;
416 0 : for (i=0;i<mod.cursor().nelements();
417 0 : i++,modIter++, delIter++, resIter++, stpIter++, masIter++) {
418 0 : if (*masIter > 0.0) {
419 0 : rHess = (*modIter + *delIter)/(1 + ggc * (*modIter + *delIter));
420 0 : gradH = -log(*modIter + *delIter) + logDef;
421 0 : gradC = -2*(*resIter);
422 0 : gradJ = gradH - alpha*gradC - beta;
423 0 : stepValue = rHess * gradJ;
424 0 : (*stpIter) = stepValue;
425 0 : GDGHH = GDGHH + gradH * rHess * gradH;
426 0 : GDGHC = GDGHC + gradH * rHess * gradC;
427 0 : GDGHF = GDGHF + gradH * rHess;
428 0 : GDGCC = GDGCC + gradC * rHess * gradC;
429 0 : GDGCF = GDGCF + gradC * rHess;
430 0 : GDGFF = GDGFF + rHess;
431 : }
432 : }
433 0 : mod.cursor().freeStorage(modStore, modDeleteIt);
434 0 : del.cursor().freeStorage(delStore, delDeleteIt);
435 0 : res.cursor().freeStorage(resStore, resDeleteIt);
436 0 : mas.cursor().freeStorage(masStore, masDeleteIt);
437 0 : stp.rwCursor().putStorage(stpStore, stpDeleteIt);
438 : }
439 : } else {
440 0 : Lattice<Float> &prior = *(cemem_ptr->itsPrior_ptr);
441 : RO_LatticeIterator<Float>
442 0 : pri( prior, LatticeStepper(prior.shape(), model.niceCursorShape()));
443 : Bool priDeleteIt;
444 0 : for (mod.reset(), del.reset(), res.reset(), pri.reset(), mas.reset(), stp.reset();
445 0 : !mod.atEnd();
446 0 : mod++, del++, res++, pri++, mas++, stp++) {
447 0 : const Float *modStore = mod.cursor().getStorage(modDeleteIt);
448 0 : const Float *delStore = del.cursor().getStorage(delDeleteIt);
449 0 : const Float *resStore = res.cursor().getStorage(resDeleteIt);
450 0 : const Float *priStore = pri.cursor().getStorage(priDeleteIt);
451 0 : const Float *masStore = mas.cursor().getStorage(masDeleteIt);
452 0 : Float *stpStore = stp.rwCursor().getStorage(stpDeleteIt);
453 0 : const Float *modIter = modStore;
454 0 : const Float *delIter = delStore;
455 0 : const Float *resIter = resStore;
456 0 : const Float *priIter = priStore;
457 0 : const Float *masIter = masStore;
458 0 : Float *stpIter = stpStore;
459 0 : for (i=0;i<mod.cursor().nelements();
460 0 : i++,modIter++, delIter++, resIter++, priIter++, stpIter++, masIter++) {
461 0 : if (*masIter > 0.0) {
462 0 : rHess = (*modIter + *delIter)/(1 + ggc* (*modIter + *delIter));
463 0 : gradH = -log( (*modIter + *delIter) / (*priIter));
464 0 : gradC = -2*(*resIter);
465 0 : gradJ = gradH - alpha*gradC - beta;
466 0 : stepValue = rHess * gradJ;
467 0 : (*stpIter) = stepValue;
468 0 : GDGHH = GDGHH + gradH * rHess * gradH;
469 0 : GDGHC = GDGHC + gradH * rHess * gradC;
470 0 : GDGHF = GDGHF + gradH * rHess;
471 0 : GDGCC = GDGCC + gradC * rHess * gradC;
472 0 : GDGCF = GDGCF + gradC * rHess;
473 0 : GDGFF = GDGFF + rHess;
474 : }
475 : }
476 0 : mod.cursor().freeStorage(modStore, modDeleteIt);
477 0 : del.cursor().freeStorage(delStore, delDeleteIt);
478 0 : res.cursor().freeStorage(resStore, resDeleteIt);
479 0 : pri.cursor().freeStorage(priStore, priDeleteIt);
480 0 : mas.cursor().freeStorage(masStore, masDeleteIt);
481 0 : stp.rwCursor().putStorage(stpStore, stpDeleteIt);
482 : }
483 0 : }
484 0 : GDG(H,H) = GDGHH;
485 0 : GDG(H,C) = GDGHC;
486 0 : GDG(H,F) = GDGHF;
487 0 : GDG(C,C) = GDGCC;
488 0 : GDG(C,F) = GDGCF;
489 0 : GDG(F,F) = GDGFF;
490 0 : GDG(H,J) = GDGHH - alpha * GDGHC - beta * GDGHF;
491 0 : GDG(C,J) = GDGHC - alpha * GDGCC - beta * GDGCF;
492 0 : GDG(F,J) = GDGHF - alpha * GDGCF - beta * GDGFF;
493 0 : GDG(J,J) = GDGHH + square(alpha) * GDGCC
494 0 : + square(beta)*GDGFF + 2*alpha*beta*GDGCF
495 0 : - 2*alpha*GDGHC - 2*beta*GDGHF;
496 :
497 :
498 0 : } else { // no mask
499 :
500 :
501 0 : Lattice<Float> &model = *(cemem_ptr->itsModel_ptr);
502 0 : Lattice<Float> &delta = *(cemem_ptr->itsDeltaModel_ptr);
503 0 : Lattice<Float> &resid = *(cemem_ptr->itsResidual_ptr);
504 0 : Lattice<Float> &step = *(cemem_ptr->itsStep_ptr);
505 :
506 0 : GDG.resize(4,4);
507 0 : GDG.set(0.0);
508 0 : Float alpha = cemem_ptr->itsAlpha;
509 0 : Float beta = cemem_ptr->itsBeta;
510 0 : Float defLev = cemem_ptr->itsDefaultLevel;
511 :
512 0 : Float ggc = 2 * alpha * cemem_ptr->itsQ;
513 : Float rHess;
514 :
515 : RO_LatticeIterator<Float>
516 0 : mod( model, LatticeStepper(model.shape(), model.niceCursorShape()));
517 : RO_LatticeIterator<Float>
518 0 : del( delta, LatticeStepper(model.shape(), model.niceCursorShape()));
519 : RO_LatticeIterator<Float>
520 0 : res( resid, LatticeStepper(resid.shape(), model.niceCursorShape()));
521 : LatticeIterator<Float>
522 0 : stp( step, LatticeStepper(step.shape(), model.niceCursorShape()));
523 : uInt i;
524 : Bool modDeleteIt;
525 : Bool delDeleteIt;
526 : Bool resDeleteIt;
527 : Bool stpDeleteIt;
528 0 : Double GDGHH = 0.0;
529 0 : Double GDGHC = 0.0;
530 0 : Double GDGHF = 0.0;
531 0 : Double GDGCC = 0.0;
532 0 : Double GDGCF = 0.0;
533 0 : Double GDGFF = 0.0;
534 0 : Double gradH = 0.0;
535 0 : Double gradC = 0.0;
536 0 : Double gradJ = 0.0;
537 : Float stepValue;
538 0 : if ( cemem_ptr->itsPrior_ptr == 0) {
539 0 : Float logDef = log(defLev);
540 0 : for (mod.reset(), del.reset(), res.reset(), stp.reset();
541 0 : !mod.atEnd();
542 0 : mod++, del++, res++, stp++) {
543 0 : const Float *modStore = mod.cursor().getStorage(modDeleteIt);
544 0 : const Float *delStore = del.cursor().getStorage(delDeleteIt);
545 0 : const Float *resStore = res.cursor().getStorage(resDeleteIt);
546 0 : Float *stpStore = stp.rwCursor().getStorage(stpDeleteIt);
547 0 : const Float *modIter = modStore;
548 0 : const Float *delIter = delStore;
549 0 : const Float *resIter = resStore;
550 0 : Float *stpIter = stpStore;
551 0 : for (i=0;i<mod.cursor().nelements();
552 0 : i++,modIter++, delIter++, resIter++, stpIter++) {
553 0 : rHess = (*modIter + *delIter)/(1 + ggc * (*modIter + *delIter));
554 0 : gradH = -log(*modIter + *delIter) + logDef;
555 0 : gradC = -2*(*resIter);
556 0 : gradJ = gradH - alpha*gradC - beta;
557 0 : stepValue = rHess * gradJ;
558 0 : (*stpIter) = stepValue;
559 0 : GDGHH = GDGHH + gradH * rHess * gradH;
560 0 : GDGHC = GDGHC + gradH * rHess * gradC;
561 0 : GDGHF = GDGHF + gradH * rHess;
562 0 : GDGCC = GDGCC + gradC * rHess * gradC;
563 0 : GDGCF = GDGCF + gradC * rHess;
564 0 : GDGFF = GDGFF + rHess;
565 : }
566 0 : mod.cursor().freeStorage(modStore, modDeleteIt);
567 0 : del.cursor().freeStorage(delStore, delDeleteIt);
568 0 : res.cursor().freeStorage(resStore, resDeleteIt);
569 0 : stp.rwCursor().putStorage(stpStore, stpDeleteIt);
570 : }
571 : } else {
572 0 : Lattice<Float> &prior = *(cemem_ptr->itsPrior_ptr);
573 : RO_LatticeIterator<Float>
574 0 : pri( prior, LatticeStepper(prior.shape(), model.niceCursorShape()));
575 : Bool priDeleteIt;
576 0 : for (mod.reset(), del.reset(), res.reset(), pri.reset(), stp.reset();
577 0 : !mod.atEnd();
578 0 : mod++, del++, res++, pri++, stp++) {
579 0 : const Float *modStore = mod.cursor().getStorage(modDeleteIt);
580 0 : const Float *delStore = del.cursor().getStorage(delDeleteIt);
581 0 : const Float *resStore = res.cursor().getStorage(resDeleteIt);
582 0 : const Float *priStore = pri.cursor().getStorage(priDeleteIt);
583 0 : Float *stpStore = stp.rwCursor().getStorage(stpDeleteIt);
584 0 : const Float *modIter = modStore;
585 0 : const Float *delIter = delStore;
586 0 : const Float *resIter = resStore;
587 0 : const Float *priIter = priStore;
588 0 : Float *stpIter = stpStore;
589 0 : for (i=0;i<mod.cursor().nelements();
590 0 : i++,modIter++, delIter++, resIter++, priIter++, stpIter++) {
591 0 : rHess = (*modIter + *delIter)/(1 + ggc* (*modIter + *delIter));
592 0 : gradH = -log( (*modIter + *delIter) / (*priIter));
593 0 : gradC = -2*(*resIter);
594 0 : gradJ = gradH - alpha*gradC - beta;
595 0 : stepValue = rHess * gradJ;
596 0 : (*stpIter) = stepValue;
597 0 : GDGHH = GDGHH + gradH * rHess * gradH;
598 0 : GDGHC = GDGHC + gradH * rHess * gradC;
599 0 : GDGHF = GDGHF + gradH * rHess;
600 0 : GDGCC = GDGCC + gradC * rHess * gradC;
601 0 : GDGCF = GDGCF + gradC * rHess;
602 0 : GDGFF = GDGFF + rHess;
603 : }
604 0 : mod.cursor().freeStorage(modStore, modDeleteIt);
605 0 : del.cursor().freeStorage(delStore, delDeleteIt);
606 0 : res.cursor().freeStorage(resStore, resDeleteIt);
607 0 : pri.cursor().freeStorage(priStore, priDeleteIt);
608 0 : stp.rwCursor().putStorage(stpStore, stpDeleteIt);
609 : }
610 0 : }
611 0 : GDG(H,H) = GDGHH;
612 0 : GDG(H,C) = GDGHC;
613 0 : GDG(H,F) = GDGHF;
614 0 : GDG(C,C) = GDGCC;
615 0 : GDG(C,F) = GDGCF;
616 0 : GDG(F,F) = GDGFF;
617 0 : GDG(H,J) = GDGHH - alpha * GDGHC - beta * GDGHF;
618 0 : GDG(C,J) = GDGHC - alpha * GDGCC - beta * GDGCF;
619 0 : GDG(F,J) = GDGHF - alpha * GDGCF - beta * GDGFF;
620 0 : GDG(J,J) = GDGHH + square(alpha) * GDGCC
621 0 : + square(beta)*GDGFF + 2*alpha*beta*GDGCF
622 0 : - 2*alpha*GDGHC - 2*beta*GDGHF;
623 :
624 0 : }
625 :
626 0 : };
627 :
628 :
629 : //----------------------------------------------------------------------
630 0 : Double IncEntropyI::formGDS()
631 : {
632 0 : Lattice<Float> &model = *(cemem_ptr->itsModel_ptr);
633 0 : Lattice<Float> &delta = *(cemem_ptr->itsDeltaModel_ptr);
634 0 : Lattice<Float> &step = *(cemem_ptr->itsStep_ptr);
635 0 : Lattice<Float> &resid = *(cemem_ptr->itsResidual_ptr);
636 :
637 0 : Float alpha = cemem_ptr->itsAlpha;
638 0 : Float beta = cemem_ptr->itsBeta;
639 0 : Float defLev = cemem_ptr->itsDefaultLevel;
640 :
641 0 : Double gds = 0;
642 :
643 0 : if (cemem_ptr->itsMask_ptr) { // Do Mask
644 0 : Lattice<Float> &mask = *(cemem_ptr->itsMask_ptr);
645 :
646 : RO_LatticeIterator<Float>
647 0 : mod( model, LatticeStepper(model.shape(), model.niceCursorShape()));
648 : RO_LatticeIterator<Float>
649 0 : del( delta, LatticeStepper(model.shape(), model.niceCursorShape()));
650 : LatticeIterator<Float>
651 0 : stp( step, LatticeStepper(model.shape(), model.niceCursorShape()));
652 : RO_LatticeIterator<Float>
653 0 : res( resid, LatticeStepper(resid.shape(), model.niceCursorShape()));
654 : RO_LatticeIterator<Float>
655 0 : mas( mask, LatticeStepper(model.shape(), model.niceCursorShape()));
656 : Bool modDeleteIt;
657 : Bool delDeleteIt;
658 : Bool masDeleteIt;
659 : Bool stpDeleteIt;
660 : Bool resDeleteIt;
661 : Bool priDeleteIt;
662 :
663 0 : if (cemem_ptr->itsPrior_ptr) { // mask AND prior
664 :
665 0 : Lattice<Float> &prior = *(cemem_ptr->itsPrior_ptr);
666 : RO_LatticeIterator<Float>
667 0 : pri( prior, LatticeStepper(model.shape(), model.niceCursorShape()));
668 :
669 :
670 0 : for (mod.reset(), del.reset(), mas.reset(), pri.reset(), res.reset(), stp.reset();
671 0 : !mod.atEnd();
672 0 : mod++, del++, mas++, pri++, res++, stp++) {
673 0 : const Float *modStore = mod.cursor().getStorage(modDeleteIt);
674 0 : const Float *delStore = del.cursor().getStorage(delDeleteIt);
675 0 : const Float *masStore = mas.cursor().getStorage(masDeleteIt);
676 0 : const Float *priStore = pri.cursor().getStorage(priDeleteIt);
677 0 : const Float *stpStore = stp.cursor().getStorage(stpDeleteIt);
678 0 : const Float *resStore = res.cursor().getStorage(resDeleteIt);
679 0 : const Float *modIter = modStore;
680 0 : const Float *delIter = delStore;
681 0 : const Float *masIter = masStore;
682 0 : const Float *priIter = priStore;
683 0 : const Float *stpIter = stpStore;
684 0 : const Float *resIter = resStore;
685 0 : for (uInt i=0;i<mod.cursor().nelements();
686 0 : i++,modIter++, delIter++, masIter++, priIter++, stpIter++,resIter++ ) {
687 0 : if (*masIter > 0.0) {
688 0 : gds += (*stpIter) * (-log( (*modIter + *delIter)/(*priIter) ) + 2.0*alpha * (*resIter) - beta);
689 : }
690 : }
691 0 : mod.cursor().freeStorage(modStore, modDeleteIt);
692 0 : mas.cursor().freeStorage(masStore, masDeleteIt);
693 0 : res.cursor().freeStorage(resStore, resDeleteIt);
694 0 : stp.cursor().freeStorage(stpStore, stpDeleteIt);
695 0 : pri.cursor().freeStorage(priStore, priDeleteIt);
696 : }
697 :
698 0 : } else { // Mask, but no prior
699 :
700 0 : Float logDefLev = log(defLev);
701 0 : for (mod.reset(), del.reset(), mas.reset(), res.reset(), stp.reset();
702 0 : !mod.atEnd();
703 0 : mod++, del++, mas++, res++, stp++) {
704 0 : const Float *modStore = mod.cursor().getStorage(modDeleteIt);
705 0 : const Float *delStore = del.cursor().getStorage(delDeleteIt);
706 0 : const Float *masStore = mas.cursor().getStorage(masDeleteIt);
707 0 : const Float *stpStore = stp.cursor().getStorage(stpDeleteIt);
708 0 : const Float *resStore = res.cursor().getStorage(resDeleteIt);
709 0 : const Float *modIter = modStore;
710 0 : const Float *delIter = delStore;
711 0 : const Float *masIter = masStore;
712 0 : const Float *stpIter = stpStore;
713 0 : const Float *resIter = resStore;
714 0 : for (uInt i=0;i<mod.cursor().nelements();i++,modIter++, delIter++, masIter++, stpIter++, resIter++) {
715 0 : if (*masIter > 0.0) {
716 0 : gds += (*stpIter) * (-log( *modIter + *delIter) +logDefLev + 2.0*alpha * (*resIter) - beta);
717 : }
718 : }
719 0 : mod.cursor().freeStorage(modStore, modDeleteIt);
720 0 : del.cursor().freeStorage(delStore, delDeleteIt);
721 0 : mas.cursor().freeStorage(masStore, masDeleteIt);
722 0 : res.cursor().freeStorage(resStore, resDeleteIt);
723 0 : stp.cursor().freeStorage(stpStore, stpDeleteIt);
724 : }
725 : }
726 :
727 0 : } else { // No Mask
728 :
729 0 : if (cemem_ptr->itsPrior_ptr == 0) {
730 0 : Float logDefLev = log(defLev);
731 0 : gds = sum( step *
732 0 : ( -log(model + delta)+logDefLev+2.0*alpha * resid -beta)).getDouble();
733 : } else {
734 0 : Lattice<Float> &prior = *(cemem_ptr->itsPrior_ptr);
735 0 : gds = sum( step *
736 0 : ( -log((model+delta)/prior) + 2.0*alpha * resid - beta)).getDouble();
737 : }
738 : }
739 :
740 0 : return gds;
741 :
742 : };
743 :
744 :
745 : //----------------------------------------------------------------------
746 0 : void IncEntropyI::infoBanner(){
747 0 : cemem_ptr->itsLog <<
748 0 : " I Entropy Flux Fit MaxResid Gradient " << LogIO::POST;
749 0 : };
750 :
751 : //----------------------------------------------------------------------
752 0 : void IncEntropyI::infoPerIteration(uInt iteration){
753 :
754 0 : Lattice<Float> &resid = *(cemem_ptr->itsResidual_ptr);
755 0 : Float maxResid = (max(abs(resid))).getFloat();
756 :
757 0 : cemem_ptr->itsLog <<
758 : (iteration +1) << " " <<
759 0 : cemem_ptr->itsEntropy << " " <<
760 0 : cemem_ptr->itsFlux << " " <<
761 0 : cemem_ptr->itsFit << " " <<
762 : maxResid << " " <<
763 0 : cemem_ptr->itsNormGrad << " " <<
764 0 : LogIO::POST;
765 :
766 0 : if (cemem_ptr->itsProgressPtr) {
767 0 : Lattice<Float> &model = *(cemem_ptr->itsModel_ptr);
768 : // Note: posMaximumResidual is not valid! ( IPosition(4,0,0,0,0) )
769 0 : cemem_ptr->itsProgressPtr->
770 0 : info(false, iteration, cemem_ptr->itsNumberIterations, model, resid,
771 0 : maxResid, IPosition(4,0,0,0,0),
772 0 : cemem_ptr->itsFlux, cemem_ptr->itsFit, cemem_ptr->itsNormGrad,
773 0 : cemem_ptr->itsEntropy);
774 : }
775 0 : };
776 :
777 :
778 : //----------------------------------------------------------------------
779 0 : Float IncEntropyI::relaxMin()
780 : {
781 0 : Float requiredModelMin = 0.1 * cemem_ptr->itsModelMin + 1.0e-9;
782 0 : return requiredModelMin;
783 : };
784 :
785 :
786 : //----------------------------------------------------------------------
787 0 : Bool IncEntropyI::testConvergence()
788 : {
789 : Bool converged;
790 0 : converged = ( (cemem_ptr->itsFit < (1.0+(cemem_ptr->itsTolerance))) &&
791 0 : (cemem_ptr->itsIteration != (cemem_ptr->itsFirstIteration)) &&
792 0 : (cemem_ptr->itsNormGrad < (cemem_ptr->itsTolerance)) );
793 0 : if (cemem_ptr->itsUseFluxConstraint) {
794 0 : converged = (converged &&
795 0 : (abs(cemem_ptr->itsFlux - cemem_ptr->itsTargetFlux)
796 0 : < cemem_ptr->itsTolerance*cemem_ptr->itsTargetFlux));
797 : }
798 :
799 0 : return converged;
800 :
801 : };
802 :
803 :
804 :
805 :
806 :
807 : //----------------------------------------------------------------------
808 0 : IncEntropyEmptiness::IncEntropyEmptiness()
809 : {
810 0 : };
811 :
812 : //----------------------------------------------------------------------
813 0 : IncEntropyEmptiness::~IncEntropyEmptiness()
814 : {
815 0 : };
816 :
817 : //----------------------------------------------------------------------
818 :
819 0 : Float IncEntropyEmptiness::formEntropy()
820 : {
821 0 : Lattice<Float> &model = *(cemem_ptr->itsModel_ptr);
822 0 : Lattice<Float> &delta = *(cemem_ptr->itsDeltaModel_ptr);
823 0 : Float flux = 0.0;
824 0 : flux = sum(model).getFloat();
825 0 : cemem_ptr->itsFlux = flux;
826 0 : LatticeExprNode myEntropyLEN;
827 0 : Float aFit = cemem_ptr->itsAFit;
828 0 : Float myEntropy = 0;
829 :
830 0 : if (cemem_ptr->itsMask_ptr) {
831 0 : Lattice<Float> &mask = *(cemem_ptr->itsMask_ptr);
832 :
833 0 : if (cemem_ptr->itsPrior_ptr == 0 ) {
834 0 : Float defLev = cemem_ptr->itsDefaultLevel;
835 0 : myEntropyLEN = - sum( mask * (log(cosh(((model + delta) - defLev)/aFit))) ) ;
836 0 : myEntropy = -aFit * myEntropyLEN.getFloat();
837 : } else {
838 0 : Lattice<Float> &prior = *(cemem_ptr->itsPrior_ptr);
839 0 : myEntropyLEN = - sum( mask * (log(cosh(((model + delta) - prior)/aFit))) ) ;
840 0 : myEntropy = -aFit * myEntropyLEN.getFloat();
841 : }
842 :
843 : } else { // No Mask
844 :
845 0 : if (cemem_ptr->itsPrior_ptr == 0 ) {
846 0 : Float defLev = cemem_ptr->itsDefaultLevel;
847 0 : myEntropyLEN = - sum( log(cosh(((model + delta) - defLev)/aFit)) ) ;
848 0 : myEntropy = -aFit * myEntropyLEN.getFloat();
849 : } else {
850 0 : Lattice<Float> &prior = *(cemem_ptr->itsPrior_ptr);
851 0 : myEntropyLEN = - sum( log(cosh(((model + delta) - prior)/aFit)) ) ;
852 0 : myEntropy = -aFit * myEntropyLEN.getFloat();
853 : }
854 : }
855 0 : return myEntropy;
856 0 : };
857 :
858 : //----------------------------------------------------------------------
859 0 : void IncEntropyEmptiness::formGDG(Matrix<Double>& GDG)
860 : {
861 0 : Lattice<Float> &delta = *(cemem_ptr->itsDeltaModel_ptr);
862 0 : Lattice<Float> &model = *(cemem_ptr->itsModel_ptr);
863 0 : Lattice<Float> &resid = *(cemem_ptr->itsResidual_ptr);
864 :
865 0 : GDG.resize(4,4);
866 0 : GDG.set(0.0);
867 0 : Float aFit = cemem_ptr->itsAFit;
868 0 : Float alpha = cemem_ptr->itsAlpha;
869 0 : Float beta = cemem_ptr->itsBeta;
870 :
871 0 : Float ggc = 2 * alpha * cemem_ptr->itsQ;
872 : Float rHess;
873 :
874 : RO_LatticeIterator<Float>
875 0 : mod( model, LatticeStepper(model.shape(), model.niceCursorShape()));
876 : RO_LatticeIterator<Float>
877 0 : del( delta, LatticeStepper(model.shape(), model.niceCursorShape()));
878 : RO_LatticeIterator<Float>
879 0 : res( resid, LatticeStepper(model.shape(), model.niceCursorShape()));
880 : uInt i;
881 : Bool modDeleteIt;
882 : Bool delDeleteIt;
883 : Bool resDeleteIt;
884 : Bool masDeleteIt;
885 0 : Double GDGHH = 0.0;
886 0 : Double GDGHC = 0.0;
887 0 : Double GDGHF = 0.0;
888 0 : Double GDGCC = 0.0;
889 0 : Double GDGCF = 0.0;
890 0 : Double GDGFF = 0.0;
891 0 : Double gradH = 0.0;
892 0 : Double gradC = 0.0;
893 :
894 0 : if (cemem_ptr->itsMask_ptr) { // Do MASK
895 :
896 0 : Lattice<Float> &mask = *(cemem_ptr->itsMask_ptr);
897 : RO_LatticeIterator<Float>
898 0 : mas( mask, LatticeStepper(model.shape(), model.niceCursorShape()));
899 :
900 0 : if (cemem_ptr->itsPrior_ptr == 0) {
901 0 : Float defLev = cemem_ptr->itsDefaultLevel;
902 0 : for (mod.reset(), del.reset(), res.reset();
903 0 : !mod.atEnd();
904 0 : mod++, del++, res++) {
905 0 : const Float *modStore = mod.cursor().getStorage(modDeleteIt);
906 0 : const Float *delStore = del.cursor().getStorage(delDeleteIt);
907 0 : const Float *resStore = res.cursor().getStorage(resDeleteIt);
908 0 : const Float *masStore = mas.cursor().getStorage(masDeleteIt);
909 0 : const Float *modIter = modStore;
910 0 : const Float *delIter = delStore;
911 0 : const Float *resIter = resStore;
912 0 : const Float *masIter = masStore;
913 0 : for (i=0;i<mod.cursor().nelements();
914 0 : i++,modIter++, delIter++, resIter++, masIter++) {
915 0 : if (*masIter > 0.0) {
916 0 : gradH = -tanh( (*modIter + *delIter - defLev)/aFit );
917 0 : rHess = 1.0/( square(1.0-gradH) /aFit + ggc) ;
918 0 : gradC = -2*(*resIter);
919 0 : GDGHH = GDGHH + gradH * rHess * gradH;
920 0 : GDGHC = GDGHC + gradH * rHess * gradC;
921 0 : GDGHF = GDGHF + gradH * rHess;
922 0 : GDGCC = GDGCC + gradC * rHess * gradC;
923 0 : GDGCF = GDGCF + gradC * rHess;
924 0 : GDGFF = GDGFF + rHess;
925 : }
926 : }
927 0 : mod.cursor().freeStorage(modStore, modDeleteIt);
928 0 : del.cursor().freeStorage(delStore, delDeleteIt);
929 0 : res.cursor().freeStorage(resStore, resDeleteIt);
930 0 : mas.cursor().freeStorage(masStore, masDeleteIt);
931 : }
932 : } else {
933 0 : Lattice<Float> &prior = *(cemem_ptr->itsPrior_ptr);
934 : RO_LatticeIterator<Float>
935 0 : pri( prior, LatticeStepper(prior.shape(), model.niceCursorShape()));
936 : Bool priDeleteIt;
937 0 : for (mod.reset(), del.reset(), res.reset(), pri.reset();
938 0 : !mod.atEnd();
939 0 : mod++, del++, res++, pri++) {
940 0 : const Float *modStore = mod.cursor().getStorage(modDeleteIt);
941 0 : const Float *delStore = del.cursor().getStorage(delDeleteIt);
942 0 : const Float *resStore = res.cursor().getStorage(resDeleteIt);
943 0 : const Float *masStore = mas.cursor().getStorage(masDeleteIt);
944 0 : const Float *priStore = pri.cursor().getStorage(priDeleteIt);
945 0 : const Float *modIter = modStore;
946 0 : const Float *delIter = delStore;
947 0 : const Float *resIter = resStore;
948 0 : const Float *masIter = masStore;
949 0 : const Float *priIter = priStore;
950 0 : for (i=0;i<mod.cursor().nelements();
951 0 : i++,modIter++, delIter++, resIter++, masIter++, priIter++) {
952 0 : if (*masIter > 0.0) {
953 0 : gradH = -tanh( (*modIter + *delIter - *priIter)/aFit );
954 0 : rHess = 1.0/( square(1.0-gradH)/aFit + ggc) ;
955 0 : gradC = -2*(*resIter);
956 0 : GDGHH = GDGHH + gradH * rHess * gradH;
957 0 : GDGHC = GDGHC + gradH * rHess * gradC;
958 0 : GDGHF = GDGHF + gradH * rHess;
959 0 : GDGCC = GDGCC + gradC * rHess * gradC;
960 0 : GDGCF = GDGCF + gradC * rHess;
961 0 : GDGFF = GDGFF + rHess;
962 : }
963 : }
964 0 : mod.cursor().freeStorage(modStore, modDeleteIt);
965 0 : del.cursor().freeStorage(delStore, delDeleteIt);
966 0 : res.cursor().freeStorage(resStore, resDeleteIt);
967 0 : mas.cursor().freeStorage(masStore, masDeleteIt);
968 0 : pri.cursor().freeStorage(priStore, priDeleteIt);
969 : }
970 0 : }
971 :
972 0 : } else { // No Mask
973 :
974 0 : if (cemem_ptr->itsPrior_ptr == 0) {
975 0 : Float defLev = cemem_ptr->itsDefaultLevel;
976 0 : for (mod.reset(), del.reset(), res.reset();
977 0 : !mod.atEnd();
978 0 : mod++, del++, res++) {
979 0 : const Float *modStore = mod.cursor().getStorage(modDeleteIt);
980 0 : const Float *delStore = del.cursor().getStorage(delDeleteIt);
981 0 : const Float *resStore = res.cursor().getStorage(resDeleteIt);
982 0 : const Float *modIter = modStore;
983 0 : const Float *delIter = delStore;
984 0 : const Float *resIter = resStore;
985 0 : for (i=0;i<mod.cursor().nelements();
986 0 : i++,modIter++, delIter++, resIter++) {
987 0 : gradH = -tanh( (*modIter + *delIter - defLev)/aFit );
988 0 : rHess = 1.0/( square(1.0-gradH) /aFit + ggc) ;
989 0 : gradC = -2*(*resIter);
990 0 : GDGHH = GDGHH + gradH * rHess * gradH;
991 0 : GDGHC = GDGHC + gradH * rHess * gradC;
992 0 : GDGHF = GDGHF + gradH * rHess;
993 0 : GDGCC = GDGCC + gradC * rHess * gradC;
994 0 : GDGCF = GDGCF + gradC * rHess;
995 0 : GDGFF = GDGFF + rHess;
996 : }
997 0 : mod.cursor().freeStorage(modStore, modDeleteIt);
998 0 : del.cursor().freeStorage(delStore, delDeleteIt);
999 0 : res.cursor().freeStorage(resStore, resDeleteIt);
1000 : }
1001 : } else {
1002 0 : Lattice<Float> &prior = *(cemem_ptr->itsPrior_ptr);
1003 : RO_LatticeIterator<Float>
1004 0 : pri( prior, LatticeStepper(prior.shape(), model.niceCursorShape()));
1005 : Bool priDeleteIt;
1006 0 : for (mod.reset(), del.reset(), res.reset(), pri.reset();
1007 0 : !mod.atEnd();
1008 0 : mod++, del++, res++, pri++) {
1009 0 : const Float *modStore = mod.cursor().getStorage(modDeleteIt);
1010 0 : const Float *delStore = del.cursor().getStorage(delDeleteIt);
1011 0 : const Float *resStore = res.cursor().getStorage(resDeleteIt);
1012 0 : const Float *priStore = pri.cursor().getStorage(priDeleteIt);
1013 0 : const Float *modIter = modStore;
1014 0 : const Float *delIter = delStore;
1015 0 : const Float *resIter = resStore;
1016 0 : const Float *priIter = priStore;
1017 0 : for (i=0;i<mod.cursor().nelements();
1018 0 : i++, modIter++, delIter++, resIter++, priIter++) {
1019 0 : gradH = -tanh( (*modIter + *delIter - *priIter)/aFit );
1020 0 : rHess = 1.0/( square(1.0-gradH)/aFit + ggc) ;
1021 0 : gradC = -2*(*resIter);
1022 0 : GDGHH = GDGHH + gradH * rHess * gradH;
1023 0 : GDGHC = GDGHC + gradH * rHess * gradC;
1024 0 : GDGHF = GDGHF + gradH * rHess;
1025 0 : GDGCC = GDGCC + gradC * rHess * gradC;
1026 0 : GDGCF = GDGCF + gradC * rHess;
1027 0 : GDGFF = GDGFF + rHess;
1028 : }
1029 0 : mod.cursor().freeStorage(modStore, modDeleteIt);
1030 0 : del.cursor().freeStorage(delStore, delDeleteIt);
1031 0 : res.cursor().freeStorage(resStore, resDeleteIt);
1032 0 : pri.cursor().freeStorage(priStore, priDeleteIt);
1033 : }
1034 0 : }
1035 :
1036 : }
1037 0 : GDG(H,H) = GDGHH;
1038 0 : GDG(H,C) = GDGHC;
1039 0 : GDG(H,F) = GDGHF;
1040 0 : GDG(C,C) = GDGCC;
1041 0 : GDG(C,F) = GDGCF;
1042 0 : GDG(F,F) = GDGFF;
1043 0 : GDG(H,J) = GDGHH - alpha * GDGHC - beta * GDGHF;
1044 0 : GDG(C,J) = GDGHC - alpha * GDGCC - beta * GDGCF;
1045 0 : GDG(F,J) = GDGHF - alpha * GDGCF - beta * GDGFF;
1046 0 : GDG(J,J) = GDGHH + square(alpha) * GDGCC
1047 0 : + square(beta)*GDGFF + 2*alpha*beta*GDGCF
1048 0 : - 2*alpha*GDGHC - 2*beta*GDGHF;
1049 0 : };
1050 :
1051 :
1052 : //----------------------------------------------------------------------
1053 0 : void IncEntropyEmptiness::formGDGStep(Matrix<Double>& GDG)
1054 : {
1055 0 : Lattice<Float> &model = *(cemem_ptr->itsModel_ptr);
1056 0 : Lattice<Float> &delta = *(cemem_ptr->itsDeltaModel_ptr);
1057 0 : Lattice<Float> &resid = *(cemem_ptr->itsResidual_ptr);
1058 0 : Lattice<Float> &step = *(cemem_ptr->itsStep_ptr);
1059 :
1060 0 : GDG.resize(4,4);
1061 0 : GDG.set(0.0);
1062 0 : Float alpha = cemem_ptr->itsAlpha;
1063 0 : Float beta = cemem_ptr->itsBeta;
1064 0 : Float aFit = cemem_ptr->itsAFit;
1065 :
1066 : // should not use Lattice Expression Language for efficiency
1067 : // using it right now to get up and running
1068 0 : Float ggc = 2 * alpha * cemem_ptr->itsQ;
1069 : Float rHess;
1070 :
1071 : RO_LatticeIterator<Float>
1072 0 : mod( model, LatticeStepper(model.shape(), model.niceCursorShape()));
1073 : RO_LatticeIterator<Float>
1074 0 : del( delta, LatticeStepper(model.shape(), model.niceCursorShape()));
1075 : RO_LatticeIterator<Float>
1076 0 : res( resid, LatticeStepper(model.shape(), model.niceCursorShape()));
1077 : LatticeIterator<Float>
1078 0 : stp( step, LatticeStepper(model.shape(), model.niceCursorShape()));
1079 : uInt i;
1080 : Bool modDeleteIt;
1081 : Bool delDeleteIt;
1082 : Bool resDeleteIt;
1083 : Bool stpDeleteIt;
1084 : Bool masDeleteIt;
1085 0 : Double GDGHH = 0.0;
1086 0 : Double GDGHC = 0.0;
1087 0 : Double GDGHF = 0.0;
1088 0 : Double GDGCC = 0.0;
1089 0 : Double GDGCF = 0.0;
1090 0 : Double GDGFF = 0.0;
1091 0 : Double gradH = 0.0;
1092 0 : Double gradC = 0.0;
1093 0 : Double gradJ = 0.0;
1094 : Float stepValue;
1095 :
1096 0 : if (cemem_ptr->itsMask_ptr) { // Mask
1097 0 : Lattice<Float> &mask = *(cemem_ptr->itsMask_ptr);
1098 : RO_LatticeIterator<Float>
1099 0 : mas( mask, LatticeStepper(model.shape(), model.niceCursorShape()));
1100 :
1101 0 : if ( cemem_ptr->itsPrior_ptr == 0) {
1102 0 : Float defLev = cemem_ptr->itsDefaultLevel;
1103 0 : for (mod.reset(), del.reset(), res.reset(), stp.reset(), mas.reset();
1104 0 : !mod.atEnd(); mod++, del++, res++, stp++, mas++) {
1105 0 : const Float *modStore = mod.cursor().getStorage(modDeleteIt);
1106 0 : const Float *delStore = del.cursor().getStorage(delDeleteIt);
1107 0 : const Float *resStore = res.cursor().getStorage(resDeleteIt);
1108 0 : const Float *masStore = mas.cursor().getStorage(masDeleteIt);
1109 0 : Float *stpStore = stp.rwCursor().getStorage(stpDeleteIt);
1110 0 : const Float *modIter = modStore;
1111 0 : const Float *delIter = delStore;
1112 0 : const Float *resIter = resStore;
1113 0 : const Float *masIter = masStore;
1114 0 : Float *stpIter = stpStore;
1115 0 : for (i=0;i<mod.cursor().nelements();
1116 0 : i++,modIter++, delIter++, resIter++, stpIter++, masIter++) {
1117 0 : if (*masIter > 0.0) {
1118 0 : gradH = -tanh( (*modIter + *delIter - defLev)/aFit );
1119 0 : rHess = 1.0/( square(1.0-gradH)/aFit + ggc) ;
1120 0 : gradC = -2*(*resIter);
1121 0 : gradJ = gradH - alpha*gradC - beta;
1122 0 : stepValue = rHess * gradJ;
1123 0 : (*stpIter) = stepValue;
1124 0 : GDGHH = GDGHH + gradH * rHess * gradH;
1125 0 : GDGHC = GDGHC + gradH * rHess * gradC;
1126 0 : GDGHF = GDGHF + gradH * rHess;
1127 0 : GDGCC = GDGCC + gradC * rHess * gradC;
1128 0 : GDGCF = GDGCF + gradC * rHess;
1129 0 : GDGFF = GDGFF + rHess;
1130 : }
1131 : }
1132 0 : mod.cursor().freeStorage(modStore, modDeleteIt);
1133 0 : del.cursor().freeStorage(delStore, delDeleteIt);
1134 0 : res.cursor().freeStorage(resStore, resDeleteIt);
1135 0 : mas.cursor().freeStorage(masStore, masDeleteIt);
1136 0 : stp.rwCursor().putStorage(stpStore, stpDeleteIt);
1137 : }
1138 : } else {
1139 0 : Lattice<Float> &prior = *(cemem_ptr->itsPrior_ptr);
1140 : RO_LatticeIterator<Float>
1141 0 : pri( prior, LatticeStepper(prior.shape(), model.niceCursorShape()));
1142 : Bool priDeleteIt;
1143 0 : for (mod.reset(), del.reset(), res.reset(), mas.reset(), pri.reset(), stp.reset();
1144 0 : !mod.atEnd(); mod++, del++, res++, mas++, pri++, stp++) {
1145 0 : const Float *modStore = mod.cursor().getStorage(modDeleteIt);
1146 0 : const Float *delStore = del.cursor().getStorage(delDeleteIt);
1147 0 : const Float *resStore = res.cursor().getStorage(resDeleteIt);
1148 0 : const Float *masStore = mas.cursor().getStorage(masDeleteIt);
1149 0 : const Float *priStore = pri.cursor().getStorage(priDeleteIt);
1150 0 : Float *stpStore = stp.rwCursor().getStorage(stpDeleteIt);
1151 0 : const Float *modIter = modStore;
1152 0 : const Float *delIter = delStore;
1153 0 : const Float *resIter = resStore;
1154 0 : const Float *masIter = masStore;
1155 0 : const Float *priIter = priStore;
1156 0 : Float *stpIter = stpStore;
1157 0 : for (i=0;i<mod.cursor().nelements();
1158 0 : i++,modIter++, delIter++, resIter++, masIter++, priIter++, stpIter++) {
1159 0 : if (*masIter > 0.0) {
1160 0 : gradH = -tanh( (*modIter + *delIter - *priIter)/aFit );
1161 0 : rHess = 1.0/( square(1.0-gradH)/aFit + ggc) ;
1162 0 : gradC = -2*(*resIter);
1163 0 : gradJ = gradH - alpha*gradC - beta;
1164 0 : stepValue = rHess * gradJ;
1165 0 : (*stpIter) = stepValue;
1166 0 : GDGHH = GDGHH + gradH * rHess * gradH;
1167 0 : GDGHC = GDGHC + gradH * rHess * gradC;
1168 0 : GDGHF = GDGHF + gradH * rHess;
1169 0 : GDGCC = GDGCC + gradC * rHess * gradC;
1170 0 : GDGCF = GDGCF + gradC * rHess;
1171 0 : GDGFF = GDGFF + rHess;
1172 : }
1173 : }
1174 0 : mod.cursor().freeStorage(modStore, modDeleteIt);
1175 0 : del.cursor().freeStorage(delStore, delDeleteIt);
1176 0 : res.cursor().freeStorage(resStore, resDeleteIt);
1177 0 : mas.cursor().freeStorage(masStore, masDeleteIt);
1178 0 : pri.cursor().freeStorage(priStore, priDeleteIt);
1179 0 : stp.rwCursor().putStorage(stpStore, stpDeleteIt);
1180 : }
1181 0 : }
1182 :
1183 0 : } else { // No Mask
1184 :
1185 0 : if ( cemem_ptr->itsPrior_ptr == 0) {
1186 0 : Float defLev = cemem_ptr->itsDefaultLevel;
1187 0 : for (mod.reset(), del.reset(), res.reset(), stp.reset();
1188 0 : !mod.atEnd(); mod++, del++, res++, stp++) {
1189 0 : const Float *modStore = mod.cursor().getStorage(modDeleteIt);
1190 0 : const Float *delStore = del.cursor().getStorage(delDeleteIt);
1191 0 : const Float *resStore = res.cursor().getStorage(resDeleteIt);
1192 0 : Float *stpStore = stp.rwCursor().getStorage(stpDeleteIt);
1193 0 : const Float *modIter = modStore;
1194 0 : const Float *delIter = delStore;
1195 0 : const Float *resIter = resStore;
1196 0 : Float *stpIter = stpStore;
1197 0 : for (i=0;i<mod.cursor().nelements();i++,modIter++, delIter++, resIter++, stpIter++) {
1198 0 : gradH = -tanh( (*modIter + *delIter - defLev)/aFit );
1199 0 : rHess = 1.0/( square(1.0-gradH)/aFit + ggc) ;
1200 0 : gradC = -2*(*resIter);
1201 0 : gradJ = gradH - alpha*gradC - beta;
1202 0 : stepValue = rHess * gradJ;
1203 0 : (*stpIter) = stepValue;
1204 0 : GDGHH = GDGHH + gradH * rHess * gradH;
1205 0 : GDGHC = GDGHC + gradH * rHess * gradC;
1206 0 : GDGHF = GDGHF + gradH * rHess;
1207 0 : GDGCC = GDGCC + gradC * rHess * gradC;
1208 0 : GDGCF = GDGCF + gradC * rHess;
1209 0 : GDGFF = GDGFF + rHess;
1210 : }
1211 0 : mod.cursor().freeStorage(modStore, modDeleteIt);
1212 0 : del.cursor().freeStorage(delStore, delDeleteIt);
1213 0 : res.cursor().freeStorage(resStore, resDeleteIt);
1214 0 : stp.rwCursor().putStorage(stpStore, stpDeleteIt);
1215 : }
1216 : } else {
1217 0 : Lattice<Float> &prior = *(cemem_ptr->itsPrior_ptr);
1218 : RO_LatticeIterator<Float>
1219 0 : pri( prior, LatticeStepper(prior.shape(), model.niceCursorShape()));
1220 : Bool priDeleteIt;
1221 0 : for (mod.reset(), del.reset(), res.reset(), pri.reset(), stp.reset();
1222 0 : !mod.atEnd(); mod++, del++, res++, pri++, stp++) {
1223 0 : const Float *modStore = mod.cursor().getStorage(modDeleteIt);
1224 0 : const Float *delStore = del.cursor().getStorage(delDeleteIt);
1225 0 : const Float *resStore = res.cursor().getStorage(resDeleteIt);
1226 0 : const Float *priStore = pri.cursor().getStorage(priDeleteIt);
1227 0 : Float *stpStore = stp.rwCursor().getStorage(stpDeleteIt);
1228 0 : const Float *modIter = modStore;
1229 0 : const Float *delIter = delStore;
1230 0 : const Float *resIter = resStore;
1231 0 : const Float *priIter = priStore;
1232 0 : Float *stpIter = stpStore;
1233 0 : for (i=0;i<mod.cursor().nelements();
1234 0 : i++,modIter++, delIter++, resIter++, priIter++, stpIter++) {
1235 0 : gradH = -tanh( (*modIter + *delIter - *priIter)/aFit );
1236 0 : rHess = 1.0/( square(1.0-gradH)/aFit + ggc) ;
1237 0 : gradC = -2*(*resIter);
1238 0 : gradJ = gradH - alpha*gradC - beta;
1239 0 : stepValue = rHess * gradJ;
1240 0 : (*stpIter) = stepValue;
1241 0 : GDGHH = GDGHH + gradH * rHess * gradH;
1242 0 : GDGHC = GDGHC + gradH * rHess * gradC;
1243 0 : GDGHF = GDGHF + gradH * rHess;
1244 0 : GDGCC = GDGCC + gradC * rHess * gradC;
1245 0 : GDGCF = GDGCF + gradC * rHess;
1246 0 : GDGFF = GDGFF + rHess;
1247 : }
1248 0 : mod.cursor().freeStorage(modStore, modDeleteIt);
1249 0 : del.cursor().freeStorage(delStore, delDeleteIt);
1250 0 : res.cursor().freeStorage(resStore, resDeleteIt);
1251 0 : pri.cursor().freeStorage(priStore, priDeleteIt);
1252 0 : stp.rwCursor().putStorage(stpStore, stpDeleteIt);
1253 : }
1254 0 : }
1255 : }
1256 0 : GDG(H,H) = GDGHH;
1257 0 : GDG(H,C) = GDGHC;
1258 0 : GDG(H,F) = GDGHF;
1259 0 : GDG(C,C) = GDGCC;
1260 0 : GDG(C,F) = GDGCF;
1261 0 : GDG(F,F) = GDGFF;
1262 0 : GDG(H,J) = GDGHH - alpha * GDGHC - beta * GDGHF;
1263 0 : GDG(C,J) = GDGHC - alpha * GDGCC - beta * GDGCF;
1264 0 : GDG(F,J) = GDGHF - alpha * GDGCF - beta * GDGFF;
1265 0 : GDG(J,J) = GDGHH + square(alpha) * GDGCC
1266 0 : + square(beta)*GDGFF + 2*alpha*beta*GDGCF
1267 0 : - 2*alpha*GDGHC - 2*beta*GDGHF;
1268 0 : };
1269 :
1270 : //----------------------------------------------------------------------
1271 0 : Double IncEntropyEmptiness::formGDS()
1272 : {
1273 0 : Lattice<Float> &model = *(cemem_ptr->itsModel_ptr);
1274 0 : Lattice<Float> &delta = *(cemem_ptr->itsDeltaModel_ptr);
1275 0 : Lattice<Float> &step = *(cemem_ptr->itsStep_ptr);
1276 0 : Lattice<Float> &resid = *(cemem_ptr->itsResidual_ptr);
1277 0 : Float alpha = cemem_ptr->itsAlpha;
1278 0 : Float beta = cemem_ptr->itsBeta;
1279 0 : Float aFit = cemem_ptr->itsAFit;
1280 :
1281 0 : Double gds = 0;
1282 :
1283 0 : if (cemem_ptr->itsMask_ptr) {
1284 0 : Lattice<Float> &mask = *(cemem_ptr->itsMask_ptr);
1285 :
1286 0 : if (cemem_ptr->itsPrior_ptr == 0) {
1287 0 : Float defLev = cemem_ptr->itsDefaultLevel;
1288 0 : gds = sum( mask * (step *
1289 0 : (-tanh((model + delta - defLev)/aFit) + 2.0*alpha * resid -beta)))
1290 0 : .getDouble();
1291 : } else {
1292 0 : Lattice<Float> &prior = *(cemem_ptr->itsPrior_ptr);
1293 0 : gds = sum( mask * (step *
1294 0 : (-tanh((model + delta - prior)/aFit) + 2.0*alpha * resid -beta)))
1295 0 : .getDouble();
1296 : }
1297 : } else { // no mask
1298 0 : if (cemem_ptr->itsPrior_ptr == 0) {
1299 0 : Float defLev = cemem_ptr->itsDefaultLevel;
1300 0 : gds = sum( step *
1301 0 : (-tanh((model + delta - defLev)/aFit) + 2.0*alpha * resid -beta))
1302 0 : .getDouble();
1303 : } else {
1304 0 : Lattice<Float> &prior = *(cemem_ptr->itsPrior_ptr);
1305 0 : gds = sum( step *
1306 0 : (-tanh((model + delta - prior)/aFit) + 2.0*alpha * resid -beta))
1307 0 : .getDouble();
1308 : }
1309 : }
1310 0 : return gds;
1311 : };
1312 :
1313 :
1314 : //----------------------------------------------------------------------
1315 0 : void IncEntropyEmptiness::infoBanner(){
1316 0 : cemem_ptr->itsLog <<
1317 0 : " I Flux Fit MaxResid Gradient " << LogIO::POST;
1318 0 : };
1319 :
1320 : //----------------------------------------------------------------------
1321 0 : void IncEntropyEmptiness::infoPerIteration(uInt iteration){
1322 :
1323 0 : Lattice<Float> &resid = *(cemem_ptr->itsResidual_ptr);
1324 0 : Float maxResid = (max(abs(resid))).getFloat();
1325 :
1326 0 : cemem_ptr->itsLog <<
1327 : (iteration +1) << " " <<
1328 0 : cemem_ptr->itsFlux << " " <<
1329 0 : cemem_ptr->itsFit << " " <<
1330 : maxResid << " " <<
1331 0 : cemem_ptr->itsNormGrad << " " <<
1332 0 : LogIO::POST;
1333 :
1334 0 : if (cemem_ptr->itsProgressPtr) {
1335 0 : Lattice<Float> &model = *(cemem_ptr->itsModel_ptr);
1336 : // Note: posMaximumResidual is not valid! ( IPosition(4,0,0,0,0) )
1337 0 : cemem_ptr->itsProgressPtr->
1338 0 : info(false, iteration, cemem_ptr->itsNumberIterations, model, resid,
1339 0 : maxResid, IPosition(4,0,0,0,0),
1340 0 : cemem_ptr->itsFlux, cemem_ptr->itsFit, cemem_ptr->itsNormGrad,
1341 0 : cemem_ptr->itsEntropy);
1342 : }
1343 0 : };
1344 :
1345 :
1346 : //----------------------------------------------------------------------
1347 0 : Float IncEntropyEmptiness::relaxMin()
1348 : {
1349 : // hey, Maximum Emptiness ain't GOT no minimum!
1350 0 : Float requiredModelMin = -1e+20;
1351 0 : return requiredModelMin;
1352 : };
1353 :
1354 :
1355 : //----------------------------------------------------------------------
1356 0 : Bool IncEntropyEmptiness::testConvergence()
1357 : {
1358 : Bool converged;
1359 0 : converged = ( (cemem_ptr->itsFit < cemem_ptr->itsTolerance) &&
1360 0 : (cemem_ptr->itsIteration != (cemem_ptr->itsFirstIteration)) &&
1361 0 : (cemem_ptr->itsNormGrad < (cemem_ptr->itsTolerance)) );
1362 0 : if (cemem_ptr->itsUseFluxConstraint) {
1363 0 : converged = (converged &&
1364 0 : (abs(cemem_ptr->itsFlux - cemem_ptr->itsTargetFlux)
1365 0 : < cemem_ptr->itsTolerance*cemem_ptr->itsTargetFlux));
1366 : }
1367 :
1368 0 : return converged;
1369 :
1370 : };
1371 :
1372 :
1373 :
1374 : } //# NAMESPACE CASA - END
1375 :
|
