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