Line data Source code
1 : //# ImageSkyModel.cc: Implementation of ImageSkyModel classes
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <casacore/casa/Arrays/ArrayMath.h>
29 : #include <casacore/casa/OS/HostInfo.h>
30 : #include <synthesis/MeasurementComponents/ImageSkyModel.h>
31 : #include <components/ComponentModels/SkyComponent.h>
32 : #include <components/ComponentModels/ComponentList.h>
33 : #include <casacore/images/Images/TempImage.h>
34 : #include <casacore/images/Images/SubImage.h>
35 : #include <casacore/images/Regions/ImageRegion.h>
36 : #include <casacore/casa/OS/File.h>
37 : #include <casacore/lattices/LEL/LatticeExpr.h>
38 : #include <casacore/lattices/Lattices/TiledLineStepper.h>
39 : #include <casacore/lattices/Lattices/LatticeStepper.h>
40 : #include <casacore/lattices/Lattices/LatticeIterator.h>
41 : #include <synthesis/MeasurementEquations/SkyEquation.h>
42 : #include <synthesis/TransformMachines/StokesImageUtil.h>
43 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
44 : #include <casacore/casa/Exceptions/Error.h>
45 : #include <casacore/casa/BasicSL/String.h>
46 : #include <casacore/casa/Utilities/Assert.h>
47 : #include <casacore/casa/OS/Directory.h>
48 : #include <casacore/tables/Tables/TableLock.h>
49 :
50 : #include <sstream>
51 :
52 : #include <casacore/casa/Logging/LogMessage.h>
53 : #include <casacore/casa/Logging/LogIO.h>
54 : #include <casacore/casa/Logging/LogSink.h>
55 :
56 : using namespace casacore;
57 : namespace casa { //# NAMESPACE CASA - BEGIN
58 :
59 : #define MEMFACTOR 8.0
60 : #if !(defined (AIPS_64B))
61 : # undef MEMFACTOR
62 : # define MEMFACTOR 18.0
63 : #endif
64 :
65 157 : ImageSkyModel::ImageSkyModel(const Int maxNumModels) :
66 : maxnmodels_p(maxNumModels),
67 : nmodels_p(0),
68 : componentList_p(0),
69 : pgplotter_p(0),
70 : displayProgress_p(false),
71 : cycleFactor_p(4.0),
72 : cycleSpeedup_p(-1.0),
73 : cycleMaxPsfFraction_p(0.8),
74 : donePSF_p(false),
75 : modified_p(true),
76 : dataPolRep_p(StokesImageUtil::CIRCULAR),
77 157 : workDirOnNFS_p(false), useMem_p(false), tileVol_p(1000000)
78 : {
79 :
80 157 : }
81 :
82 0 : void ImageSkyModel::setMaxNumberModels(const Int maxNumModels) {
83 0 : maxnmodels_p = maxNumModels;
84 0 : }
85 :
86 105 : Bool ImageSkyModel::add(ComponentList& compList)
87 : {
88 105 : if(componentList_p==0) {
89 105 : componentList_p=new ComponentList(compList);
90 105 : return true;
91 : }
92 0 : return false;
93 : }
94 :
95 0 : Bool ImageSkyModel::updatemodel(ComponentList& compList)
96 : {
97 0 : if(componentList_p!=0) {
98 0 : delete componentList_p;
99 0 : componentList_p=0;
100 : }
101 :
102 0 : componentList_p=new ComponentList(compList);
103 0 : modified_p=true;
104 0 : return true;
105 : }
106 :
107 :
108 0 : Bool ImageSkyModel::updatemodel(const Int thismodel, ImageInterface<Float>& image){
109 0 : if(nmodels_p < thismodel)
110 0 : throw(AipsError("Programming error " + String::toString(thismodel) + " is larger than the number of models"));
111 0 : image_p[thismodel]=ℑ
112 0 : AlwaysAssert(image_p[thismodel], AipsError);
113 0 : image_p[thismodel]->setUnits(Unit("Jy/pixel"));
114 0 : modified_p=true;
115 0 : return true;
116 : }
117 :
118 57 : Int ImageSkyModel::add(ImageInterface<Float>& image, const Int maxNumXfr)
119 : {
120 57 : Int thismodel=nmodels_p;
121 57 : nmodels_p++;
122 : /////Avoid using nfs file if memory is available
123 : try{
124 63 : Directory eldir(image.name());
125 57 : workDirOnNFS_p=eldir.isNFSMounted();
126 : //cerr << "Using NFS " << workDirOnNFS_p << endl;
127 : }
128 6 : catch(...){
129 6 : workDirOnNFS_p=false;
130 : }
131 :
132 57 : if(nmodels_p>maxnmodels_p) maxnmodels_p=nmodels_p;
133 :
134 57 : maxNumXFR_p=maxNumXfr;
135 :
136 57 : image_p.resize(nmodels_p);
137 57 : cimage_p.resize(nmodels_p);
138 57 : cxfr_p.resize(nmodels_p*maxNumXFR_p);
139 57 : residual_p.resize(nmodels_p);
140 57 : residualImage_p.resize(nmodels_p);
141 57 : gS_p.resize(nmodels_p);
142 57 : psf_p.resize(nmodels_p);
143 57 : ggS_p.resize(nmodels_p);
144 57 : fluxScale_p.resize(nmodels_p);
145 57 : work_p.resize(nmodels_p);
146 57 : deltaimage_p.resize(nmodels_p);
147 57 : solve_p.resize(nmodels_p);
148 57 : doFluxScale_p.resize(nmodels_p);
149 57 : weight_p.resize(nmodels_p);
150 57 : beam_p.resize(nmodels_p);
151 :
152 57 : image_p[thismodel]=0;
153 57 : cimage_p[thismodel]=0;
154 57 : residual_p[thismodel]=0;
155 :
156 2442 : for (Int numXFR=0;numXFR<maxNumXFR_p;numXFR++) {
157 2385 : cxfr_p[thismodel*maxNumXFR_p+numXFR]=0;
158 : }
159 57 : residualImage_p[thismodel]=0;
160 57 : gS_p[thismodel]=0;
161 57 : psf_p[thismodel]=0;
162 57 : ggS_p[thismodel]=0;
163 57 : fluxScale_p[thismodel]=0;
164 57 : work_p[thismodel]=0;
165 57 : deltaimage_p[thismodel]=0;
166 57 : solve_p[thismodel]=true;
167 57 : doFluxScale_p[thismodel]=false;
168 57 : weight_p[thismodel]=0;
169 57 : beam_p[thismodel]=0;
170 :
171 : // Initialize image
172 57 : image_p[thismodel]=ℑ
173 57 : AlwaysAssert(image_p[thismodel], AipsError);
174 57 : image_p[thismodel]->setUnits(Unit("Jy/pixel"));
175 57 : donePSF_p=false;
176 57 : return thismodel;
177 : }
178 :
179 0 : Bool ImageSkyModel::addResidual(Int thismodel, ImageInterface<Float>& residual)
180 : {
181 0 : LogIO os(LogOrigin("ImageSkyModel", "addResidual"));
182 0 : if(thismodel>=nmodels_p||thismodel<0) {
183 0 : os << LogIO::SEVERE << "Illegal model slot" << thismodel << LogIO::POST;
184 0 : return false;
185 : }
186 0 : residual_p[thismodel] = &residual;
187 0 : AlwaysAssert(residual_p[thismodel], AipsError);
188 0 : if(residualImage_p[thismodel]) delete residualImage_p[thismodel];
189 0 : residualImage_p[thismodel]=0;
190 0 : return true;
191 : }
192 :
193 0 : ImageSkyModel::ImageSkyModel(const ImageSkyModel& other) : SkyModel() {
194 0 : operator=(other);
195 0 : };
196 :
197 157 : ImageSkyModel::~ImageSkyModel() {
198 157 : if(componentList_p) delete componentList_p; componentList_p=0;
199 214 : for (Int thismodel=0;thismodel<nmodels_p;thismodel++) {
200 57 : if(residualImage_p[thismodel]) {
201 : //(residualImage_p[thismodel])->clearCache();
202 : //(residualImage_p[thismodel])->tempClose();
203 0 : delete residualImage_p[thismodel];
204 : }
205 57 : residualImage_p[thismodel]=0;
206 57 : if(cimage_p[thismodel]){
207 : //(cimage_p[thismodel])->clearCache();
208 : //(cimage_p[thismodel])->tempClose();
209 49 : delete cimage_p[thismodel];
210 : }
211 57 : cimage_p[thismodel]=0;
212 2442 : for(Int numXFR=0;numXFR<maxNumXFR_p;numXFR++) {
213 2385 : if(cxfr_p[thismodel*maxNumXFR_p+numXFR])
214 0 : delete cxfr_p[thismodel*maxNumXFR_p+numXFR];
215 2385 : cxfr_p[thismodel*maxNumXFR_p+numXFR]=0;
216 : }
217 57 : if(gS_p[thismodel]) delete gS_p[thismodel]; gS_p[thismodel]=0;
218 57 : if(psf_p[thismodel]) delete psf_p[thismodel]; psf_p[thismodel]=0;
219 57 : if(ggS_p[thismodel]) delete ggS_p[thismodel]; ggS_p[thismodel]=0;
220 57 : if(fluxScale_p[thismodel]) delete fluxScale_p[thismodel]; fluxScale_p[thismodel]=0;
221 57 : if(work_p[thismodel]) delete work_p[thismodel]; work_p[thismodel]=0;
222 57 : if(deltaimage_p[thismodel]) delete deltaimage_p[thismodel]; deltaimage_p[thismodel]=0;
223 57 : if(weight_p[thismodel]) delete weight_p[thismodel]; weight_p[thismodel]=0;
224 57 : if(beam_p[thismodel]) delete beam_p[thismodel]; beam_p[thismodel]=0;
225 : }
226 157 : };
227 :
228 0 : ImageSkyModel& ImageSkyModel::operator=(const ImageSkyModel& other) {
229 0 : if(this!=&other) {
230 0 : componentList_p=other.componentList_p;
231 0 : pgplotter_p=other.pgplotter_p;
232 0 : donePSF_p=other.donePSF_p;
233 0 : nmodels_p=other.nmodels_p;
234 0 : for (Int thismodel=0;thismodel<nmodels_p;thismodel++) {
235 0 : image_p[thismodel]=other.image_p[thismodel];
236 0 : cimage_p[thismodel]=other.cimage_p[thismodel];
237 0 : for(Int numXFR=0;numXFR<maxNumXFR_p;numXFR++) {
238 0 : cxfr_p[thismodel*maxNumXFR_p+numXFR]=
239 0 : other.cxfr_p[thismodel*maxNumXFR_p+numXFR];
240 : }
241 0 : residual_p[thismodel]=other.residual_p[thismodel];
242 0 : residualImage_p[thismodel]=other.residualImage_p[thismodel];
243 0 : gS_p[thismodel]=other.gS_p[thismodel];
244 0 : ggS_p[thismodel]=other.ggS_p[thismodel];
245 0 : fluxScale_p[thismodel]=other.fluxScale_p[thismodel];
246 0 : work_p[thismodel]=other.work_p[thismodel];
247 0 : deltaimage_p[thismodel]=other.deltaimage_p[thismodel];
248 0 : psf_p[thismodel]=other.psf_p[thismodel];
249 0 : weight_p[thismodel]=other.weight_p[thismodel];
250 0 : beam_p[thismodel]=other.beam_p[thismodel];
251 : }
252 0 : sumwt_p=other.sumwt_p;
253 0 : chisq_p=other.chisq_p;
254 : };
255 0 : return *this;
256 : }
257 :
258 : // Make the PSF image for each model
259 14 : void ImageSkyModel::makeApproxPSFs(SkyEquation& se) {
260 28 : LogIO os(LogOrigin("ImageSkyModel", "makeApproxPSFs"));
261 :
262 14 : if(!donePSF_p){
263 28 : for (Int thismodel=0;thismodel<nmodels_p;thismodel++) {
264 : //make sure the psf images are made
265 14 : PSF(thismodel);
266 : }
267 14 : se.makeApproxPSF(psf_p);
268 28 : for (Int thismodel=0;thismodel<nmodels_p;thismodel++) {
269 : //beam(thismodel) = 0.0;
270 14 : if(!StokesImageUtil::FitGaussianPSF(PSF(thismodel),
271 : beam(thismodel))) {
272 0 : os << "Beam fit failed: using default" << LogIO::POST;
273 : }
274 14 : if(nmodels_p > 1)
275 0 : os << LogIO::NORMAL << "Model " << thismodel+1 << ": "; // Loglevel INFO
276 : os << LogIO::NORMAL // Loglevel INFO
277 14 : << "bmaj: " << abs((beam(thismodel)(0,0)).getMajor("arcsec"))
278 14 : << "\", bmin: " << abs((beam(thismodel)(0,0)).getMinor("arcsec"))
279 28 : << "\", bpa: " << ((beam(thismodel)(0,0)).getPA(Unit("deg")))
280 70 : << " deg" << LogIO::POST;
281 : }
282 : }
283 14 : donePSF_p=true;
284 14 : }
285 :
286 14 : void ImageSkyModel::initializeGradients() {
287 14 : sumwt_p=0.0;
288 14 : chisq_p=0.0;
289 28 : for (Int thismodel=0;thismodel<nmodels_p;thismodel++) {
290 14 : cImage(thismodel).set(Complex(0.0));
291 14 : gS(thismodel).set(0.0);
292 14 : ggS(thismodel).set(0.0);
293 : }
294 14 : }
295 :
296 0 : Bool ImageSkyModel::makeNewtonRaphsonStep(SkyEquation& se, Bool incremental, Bool modelToMS)
297 : {
298 0 : LogIO os(LogOrigin("ImageSkyModel", "makeNewtonRaphsonStep"));
299 0 : se.gradientsChiSquared(incremental, modelToMS);
300 : // os << "Image normalization moved to CSE!!" << LogIO::WARN << LogIO::POST;
301 :
302 : // Now for each model, we find the recommended step
303 0 : LatticeExpr<Float> le;
304 0 : if(numberOfModels()>0) {
305 0 : for(Int thismodel=0;thismodel<nmodels_p;thismodel++) {
306 0 : if(isSolveable(thismodel)) {
307 : // SB
308 : // Use the following code without the CSE::tmpNormalizeImage()
309 0 : if (isImageNormalized()) le = LatticeExpr<Float>(gS(thismodel));
310 0 : else le = (iif(ggS(thismodel)>(0.0), -gS(thismodel)/ggS(thismodel), 0.0));
311 :
312 : // Use the following code when using CSE::tmpNormalizeImage()
313 : // le = LatticeExpr<Float>(gS(thismodel));
314 0 : residual(thismodel).copyData(le);
315 : }
316 : }
317 : }
318 0 : modified_p=false;
319 0 : return true;
320 : }
321 :
322 : // Simply finds residual image: i.e. Dirty Image if started with
323 : // zero'ed image. We work from corrected visibilities only!
324 0 : Bool ImageSkyModel::solve(SkyEquation& se) {
325 0 : return solveResiduals(se);
326 : }
327 :
328 : // Simply finds residual image: i.e. Dirty Image if started with
329 : // zero'ed image. We work from corrected visibilities only!
330 0 : Bool ImageSkyModel::solveResiduals(SkyEquation& se, Bool modelToMS) {
331 0 : makeNewtonRaphsonStep(se, false, modelToMS);
332 0 : return true;
333 : }
334 :
335 : // Return residual image: to be used by callers when it is known that the
336 : // current residual image is correct
337 0 : ImageInterface<Float>& ImageSkyModel::getResidual(Int model) {
338 0 : return ImageSkyModel::residual(model);
339 : }
340 :
341 66 : Bool ImageSkyModel::isEmpty(Int model)
342 : {
343 132 : const IPosition tileShape = image(model).niceCursorShape();
344 132 : TiledLineStepper ls(image(model).shape(), tileShape, 0);
345 132 : RO_LatticeIterator<Float> li(image(model), ls);
346 :
347 4057 : for(li.reset();!li.atEnd();li++) {
348 4049 : if(max(abs(li.cursor()))>0.0) return false;
349 : }
350 8 : return true;
351 : }
352 :
353 289 : ImageInterface<Float>& ImageSkyModel::image(Int model)
354 : {
355 289 : AlwaysAssert(nmodels_p>0, AipsError);
356 289 : AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
357 289 : AlwaysAssert(image_p[model], AipsError);
358 289 : return *image_p[model];
359 : };
360 :
361 :
362 :
363 351 : ImageInterface<Complex>& ImageSkyModel::cImage(Int model)
364 : {
365 351 : AlwaysAssert(nmodels_p>0, AipsError);
366 351 : AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
367 :
368 : //if(model>0&&(cimage_p[model-1])) cimage_p[model-1]->tempClose();
369 :
370 351 : Double memoryMB=HostInfo::memoryFree()/1024/(MEMFACTOR*maxnmodels_p);
371 351 : if(cimage_p[model]==0) {
372 98 : Vector<Int> whichStokes(0);
373 98 : IPosition cimageShape;
374 49 : cimageShape=image_p[model]->shape();
375 :
376 : CoordinateSystem cimageCoord =
377 : StokesImageUtil::CStokesCoord(//cimageShape,
378 98 : image_p[model]->coordinates(),
379 : whichStokes,
380 98 : dataPolRep_p);
381 : // StokesImageUtil::CIRCULAR);
382 :
383 : /* STOKESDBG */ //cout << "ImageSkyModel::CImage : Correlation Planes 'whichStokes' : " << whichStokes << endl;
384 49 : cimageShape(2)=whichStokes.nelements();
385 :
386 : // Now set up the tile size, here we guess only
387 : // IPosition tileShape(4, min(32, cimageShape(0)), min(32, cimageShape(1)),
388 : // min(4, cimageShape(2)), min(32, cimageShape(3)));
389 :
390 : //If tempImage is going to use disk based image ...try to respect the tile shape of
391 : //of original model image
392 :
393 : TempImage<Complex>* cimagePtr =
394 98 : new TempImage<Complex> (TiledShape(cimageShape, tileShape(model)),
395 : //IPosition(4, min(image_p[model]->shape()(0), 1000), min(image_p[model]->shape()(1), 1000), 1, 1)),
396 49 : cimageCoord, (workDirOnNFS_p|| useMem_p) ? memoryMB : 0);
397 49 : cimagePtr->set(0.0);
398 49 : cimagePtr->setMaximumCacheSize(min(cimagePtr->shape().product()/10, 1000000));
399 : // cimagePtr->setMaximumCacheSize(cacheSize(model));
400 49 : AlwaysAssert(cimagePtr, AipsError);
401 49 : cimage_p[model] = cimagePtr;
402 49 : cimage_p[model]->setMiscInfo(image_p[model]->miscInfo());
403 49 : cimagePtr->clearCache();
404 : }
405 351 : return *cimage_p[model];
406 : };
407 :
408 :
409 0 : Bool ImageSkyModel::hasXFR(Int model)
410 : {
411 0 : return (cxfr_p[model]);
412 : }
413 :
414 0 : ImageInterface<Complex>& ImageSkyModel::XFR(Int model, Int numXFR)
415 : {
416 :
417 0 : AlwaysAssert(nmodels_p>0, AipsError);
418 0 : AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
419 0 : if(numXFR==maxNumXFR_p){
420 0 : ++maxNumXFR_p;
421 0 : cxfr_p.resize(nmodels_p*maxNumXFR_p, true);
422 : //initialize the extra pointers to 0
423 0 : for (Int k = 0; k <nmodels_p; ++k){
424 0 : cxfr_p[nmodels_p*(maxNumXFR_p-1)+k]=0;
425 : }
426 : }
427 0 : AlwaysAssert(numXFR<maxNumXFR_p, AipsError);
428 0 : Double memoryMB=HostInfo::memoryTotal(true)/1024/(MEMFACTOR*maxnmodels_p);
429 0 : if(cxfr_p[model*maxNumXFR_p+numXFR]==0) {
430 :
431 0 : TempImage<Complex>* cxfrPtr = 0;
432 0 : if (imageRegion_p) {
433 : // use the imageRegion_p to define the image size;
434 : // this will usually be related to the primary beam size,
435 : // but is specified outside of the realm of the ImageSkyModel
436 :
437 : // cout << "ISM::XFR() Using ImageRegion to make smaller XFR" << endl;
438 :
439 0 : SubImage<Complex> sub( cImage(model), *imageRegion_p );
440 0 : cxfrPtr =
441 0 : new TempImage<Complex> (IPosition(sub.ndim(),
442 0 : sub.shape()(0),
443 0 : sub.shape()(1),
444 0 : sub.shape()(2),
445 0 : sub.shape()(3)),
446 0 : sub.coordinates(),
447 0 : workDirOnNFS_p ? memoryMB : 0);
448 : // cout << "ISM::XFR() shape = " << cxfrPtr->shape() << endl;
449 :
450 : } else {
451 : // use default (ie, full) image size
452 0 : cxfrPtr =
453 0 : new TempImage<Complex> (IPosition(cImage(model).ndim(),
454 0 : cImage(model).shape()(0),
455 0 : cImage(model).shape()(1),
456 0 : cImage(model).shape()(2),
457 0 : cImage(model).shape()(3)),
458 0 : cImage(model).coordinates(),
459 0 : workDirOnNFS_p ? memoryMB : 0);
460 : }
461 0 : AlwaysAssert(cxfrPtr, AipsError);
462 : //cxfrPtr->setMaximumCacheSize(cxfrPtr->shape().product()/2);
463 0 : cxfrPtr->setMaximumCacheSize(cacheSize(model));
464 0 : cxfrPtr->clearCache();
465 0 : cxfr_p[model*maxNumXFR_p+numXFR] = cxfrPtr;
466 : }
467 0 : AlwaysAssert(cxfr_p[model*maxNumXFR_p+numXFR], AipsError);
468 0 : return *cxfr_p[model*maxNumXFR_p+numXFR];
469 : };
470 :
471 :
472 56 : ImageInterface<Float>& ImageSkyModel::PSF(Int model)
473 : {
474 56 : AlwaysAssert(nmodels_p>0, AipsError);
475 56 : AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
476 :
477 : //if(model>0&&(psf_p[model-1])) psf_p[model-1]->tempClose();
478 :
479 56 : Double memoryMB=HostInfo::memoryTotal(true)/1024/(MEMFACTOR*maxnmodels_p);
480 56 : if(psf_p[model]==0) {
481 : TempImage<Float>* psfPtr =
482 28 : new TempImage<Float> (TiledShape(image_p[model]->shape(), tileShape(model)),
483 14 : image_p[model]->coordinates(),
484 14 : (workDirOnNFS_p || useMem_p) ? memoryMB : 0);
485 :
486 14 : AlwaysAssert(psfPtr, AipsError);
487 14 : psfPtr->set(0.0);
488 : //psfPtr->setMaximumCacheSize(psfPtr->shape().product()/10);
489 14 : psfPtr->setMaximumCacheSize(cacheSize(model));
490 14 : psfPtr->clearCache();
491 14 : psf_p[model] = psfPtr;
492 : }
493 56 : return *psf_p[model];
494 : };
495 :
496 :
497 0 : ImageInterface<Float>& ImageSkyModel::residual(Int model) {
498 0 : AlwaysAssert(nmodels_p>0, AipsError);
499 0 : AlwaysAssert(nmodels_p>0, AipsError);
500 0 : AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
501 0 : if(residual_p[model]) {
502 0 : if(residualImage_p[model]) delete residualImage_p[model];
503 0 : residualImage_p[model]=0;
504 0 : return *residual_p[model];
505 : }
506 : else {
507 0 : if(residualImage_p[model]==0) {
508 0 : Double memoryMB=HostInfo::memoryFree()/1024.0/(MEMFACTOR*maxnmodels_p);
509 :
510 : TempImage<Float>* tempImagePtr =
511 0 : new TempImage<Float> (TiledShape(image_p[model]->shape(),
512 0 : image_p[model]->niceCursorShape()),
513 : //IPosition(4, min(image_p[model]->shape()(0), 1000), min(image_p[model]->shape()(1), 1000), 1, 1)),
514 0 : image_p[model]->coordinates(), memoryMB);
515 :
516 0 : AlwaysAssert(tempImagePtr, AipsError);
517 :
518 : //tempImagePtr->setMaximumCacheSize(tempImagePtr->shape().product()/10);
519 0 : tempImagePtr->setMaximumCacheSize(cacheSize(model));
520 0 : tempImagePtr->set(0.0);
521 0 : tempImagePtr->clearCache();
522 0 : residualImage_p[model] = tempImagePtr;
523 : }
524 0 : return *residualImage_p[model];
525 : }
526 : }
527 :
528 56 : ImageInterface<Float>& ImageSkyModel::gS(Int model)
529 : {
530 56 : AlwaysAssert(nmodels_p>0, AipsError);
531 56 : AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
532 :
533 : //if(model>0&&(gS_p[model-1])) gS_p[model-1]->tempClose();
534 :
535 56 : if(gS_p[model]==0) {
536 14 : Double memoryMB=HostInfo::memoryTotal(true)/1024/(MEMFACTOR*maxnmodels_p);
537 :
538 : TempImage<Float>* gSPtr =
539 28 : new TempImage<Float> (TiledShape(image_p[model]->shape(), tileShape(model)),
540 14 : image_p[model]->coordinates(), (workDirOnNFS_p || useMem_p) ? memoryMB : 0);
541 :
542 : //gSPtr->setMaximumCacheSize(gSPtr->shape().product()/10);
543 14 : gSPtr->setMaximumCacheSize(cacheSize(model));
544 14 : gSPtr->set(0.0);
545 14 : gSPtr->clearCache();
546 14 : AlwaysAssert(gSPtr, AipsError);
547 14 : gS_p[model] = gSPtr;
548 : }
549 56 : return *gS_p[model];
550 : };
551 :
552 126 : ImageInterface<Float>& ImageSkyModel::ggS(Int model)
553 : {
554 126 : AlwaysAssert(nmodels_p>0, AipsError);
555 :
556 : //if(model>0&&(ggS_p[model-1])) ggS_p[model-1]->tempClose();
557 :
558 126 : if(ggS_p[model]==0) {
559 14 : Double memoryMB=HostInfo::memoryTotal(true)/1024/(MEMFACTOR*maxnmodels_p);
560 : TempImage<Float>* ggSPtr =
561 28 : new TempImage<Float> (TiledShape(image_p[model]->shape(), tileShape(model)),
562 14 : image_p[model]->coordinates(),
563 14 : (workDirOnNFS_p || useMem_p) ? memoryMB : 0);
564 :
565 14 : AlwaysAssert(ggSPtr, AipsError);
566 : //ggSPtr->setMaximumCacheSize(ggSPtr->shape().product()/10);
567 14 : ggSPtr->setMaximumCacheSize(cacheSize(model));
568 14 : ggSPtr->set(0.0);
569 14 : ggSPtr->clearCache();
570 14 : ggS_p[model] = ggSPtr;
571 : }
572 126 : AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
573 126 : return *ggS_p[model];
574 : };
575 :
576 0 : ImageInterface<Float>& ImageSkyModel::fluxScale(Int model)
577 : {
578 0 : AlwaysAssert(nmodels_p>0, AipsError);
579 :
580 : // if(model>0&&(fluxScale_p[model-1])) fluxScale_p[model-1]->tempClose();
581 :
582 0 : if(fluxScale_p[model]==0) {
583 0 : Double memoryMB=HostInfo::memoryTotal(true)/1024/(MEMFACTOR*maxnmodels_p);
584 :
585 : TempImage<Float>* fluxScalePtr =
586 0 : new TempImage<Float> (TiledShape(image_p[model]->shape(), tileShape(model)),
587 0 : image_p[model]->coordinates(),
588 0 : (workDirOnNFS_p || useMem_p) ? memoryMB : 0);
589 :
590 0 : AlwaysAssert(fluxScalePtr, AipsError);
591 : //fluxScalePtr->setMaximumCacheSize(fluxScalePtr->shape().product()/10);
592 0 : fluxScalePtr->setMaximumCacheSize(cacheSize(model));
593 0 : fluxScalePtr->set(0.0);
594 0 : fluxScalePtr->clearCache();
595 0 : fluxScale_p[model] = fluxScalePtr;
596 : // Set default value to avoid a nasty side effect elsewhere
597 0 : fluxScale_p[model]->set(1.0);
598 : }
599 0 : AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
600 0 : if( isImageNormalized() == false )
601 : {
602 0 : mandateFluxScale(model);
603 : }
604 0 : return *fluxScale_p[model];
605 : };
606 :
607 28 : ImageInterface<Float>& ImageSkyModel::work(Int model)
608 : {
609 28 : AlwaysAssert(nmodels_p>0, AipsError);
610 28 : AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
611 :
612 : // if(model>0&&(work_p[model-1])) work_p[model-1]->tempClose();
613 :
614 28 : if(work_p[model]==0) {
615 14 : Double memoryMB=HostInfo::memoryTotal(true)/1024/(MEMFACTOR*maxnmodels_p);
616 :
617 : TempImage<Float>* workPtr =
618 28 : new TempImage<Float> (TiledShape(image_p[model]->shape(),tileShape(model)),
619 14 : image_p[model]->coordinates(),
620 14 : (workDirOnNFS_p || useMem_p) ? memoryMB : 0);
621 :
622 14 : AlwaysAssert(workPtr, AipsError);
623 : //workPtr->setMaximumCacheSize(workPtr->shape().product()/10);
624 14 : workPtr->setMaximumCacheSize(cacheSize(model));
625 14 : workPtr->set(0.0);
626 14 : workPtr->clearCache();
627 14 : work_p[model] = workPtr;
628 : }
629 28 : return *work_p[model];
630 : };
631 :
632 2 : ImageInterface<Float>& ImageSkyModel::deltaImage(Int model)
633 : {
634 2 : AlwaysAssert(nmodels_p>0, AipsError);
635 2 : AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
636 :
637 : // if(model>0&&(deltaimage_p[model-1])) deltaimage_p[model-1]->tempClose();
638 :
639 2 : if(deltaimage_p[model]==0) {
640 1 : Double memoryMB=HostInfo::memoryTotal(true)/1024/(MEMFACTOR*maxnmodels_p);
641 : TempImage<Float>* deltaimagePtr =
642 2 : new TempImage<Float> (TiledShape(image_p[model]->shape(),tileShape(model)),
643 1 : image_p[model]->coordinates(),
644 1 : (workDirOnNFS_p || useMem_p) ? memoryMB : 0);
645 :
646 1 : AlwaysAssert(deltaimagePtr, AipsError);
647 : //deltaimagePtr->setMaximumCacheSize(deltaimagePtr->shape().product()/10);
648 1 : deltaimagePtr->setMaximumCacheSize(cacheSize(model));
649 1 : deltaimagePtr->set(0.0);
650 1 : deltaimagePtr->clearCache();
651 1 : deltaimage_p[model] = deltaimagePtr;
652 : }
653 2 : return *deltaimage_p[model];
654 : };
655 :
656 0 : Matrix<Float>& ImageSkyModel::weight(Int model)
657 : {
658 0 : AlwaysAssert(nmodels_p>0, AipsError);
659 0 : AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
660 0 : if(weight_p[model]==0) {
661 0 : weight_p[model] = new Matrix<Float>(1,1);
662 0 : AlwaysAssert(weight_p[model], AipsError);
663 0 : *weight_p[model]=Float(0.0);
664 : }
665 0 : return *weight_p[model];
666 : };
667 :
668 56 : ImageBeamSet& ImageSkyModel::beam(Int model)
669 : {
670 56 : AlwaysAssert(nmodels_p>0, AipsError);
671 56 : AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
672 56 : if(beam_p[model]==0) {
673 14 : beam_p[model] = new ImageBeamSet();
674 :
675 14 : AlwaysAssert(beam_p[model], AipsError);
676 : //*beam_p[model] = Float(-1.0);
677 : }
678 56 : return *beam_p[model];
679 : };
680 :
681 :
682 0 : Bool ImageSkyModel::free(Int model)
683 : {
684 0 : AlwaysAssert(nmodels_p>0, AipsError);
685 0 : AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
686 0 : Bool previous=solve_p[model];
687 0 : solve_p[model]=true;
688 0 : return previous;
689 : };
690 :
691 0 : Bool ImageSkyModel::fix(Int model)
692 : {
693 0 : AlwaysAssert(nmodels_p>0, AipsError);
694 0 : AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
695 0 : Bool previous=solve_p[model];
696 0 : solve_p[model]=false;
697 0 : return previous;
698 : };
699 :
700 :
701 0 : Bool ImageSkyModel::isSolveable(Int model)
702 : {
703 0 : AlwaysAssert(nmodels_p>0, AipsError);
704 0 : AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
705 0 : return solve_p[model];
706 : };
707 :
708 :
709 70 : Bool ImageSkyModel::doFluxScale(Int model)
710 : {
711 70 : AlwaysAssert(nmodels_p>0, AipsError);
712 70 : AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
713 70 : return doFluxScale_p[model];
714 : };
715 :
716 0 : void ImageSkyModel::mandateFluxScale(Int model)
717 : {
718 0 : AlwaysAssert(nmodels_p>0, AipsError);
719 0 : AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
720 0 : doFluxScale_p[model]=true;
721 0 : };
722 :
723 143 : Bool ImageSkyModel::hasComponentList()
724 : {
725 143 : return (componentList_p);
726 : }
727 :
728 5245 : ComponentList& ImageSkyModel::componentList()
729 : {
730 5245 : AlwaysAssert(componentList_p, AipsError);
731 5245 : return *componentList_p;
732 : }
733 :
734 57 : Long ImageSkyModel::cacheSize(Int model){
735 57 : AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
736 : // Long cachesize=(image_p[model]->shape()[0])*(image_p[model]->shape()[1])/4;
737 57 : Long cachesize=(image_p[model]->shape().product());
738 : // return Long((HostInfo::memoryFree()/(sizeof(Float)*16)*1024));
739 :
740 57 : return cachesize;
741 : }
742 :
743 106 : IPosition ImageSkyModel::tileShape(Int model){
744 106 : AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
745 106 : Int x=Int(sqrt(Double(tileVol_p)));
746 212 : IPosition shp(4, min(image_p[model]->shape()(0), x), min(image_p[model]->shape()(1), x), 1, 1);
747 :
748 106 : return shp;
749 : }
750 :
751 124 : void ImageSkyModel::setMemoryUse(Bool useMem){
752 124 : useMem_p=useMem;
753 124 : }
754 0 : void ImageSkyModel::setTileVol(Int tileVol){
755 0 : tileVol_p=tileVol;
756 0 : }
757 :
758 : } //# NAMESPACE CASA - END
759 :
|
