LCOV - code coverage report
Current view: top level - synthesis/MeasurementComponents - ImageSkyModel.cc (source / functions) Hit Total Coverage
Test: casacpp_coverage.info Lines: 221 428 51.6 %
Date: 2024-10-04 18:58:15 Functions: 20 42 47.6 %

          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: 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 <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          36 : ImageSkyModel::ImageSkyModel(const Int maxNumModels) :
      66          36 :   maxnmodels_p(maxNumModels), 
      67          36 :   nmodels_p(0), 
      68          36 :   componentList_p(0), 
      69          36 :   pgplotter_p(0),
      70          36 :   displayProgress_p(false),
      71          36 :   cycleFactor_p(4.0),
      72          36 :   cycleSpeedup_p(-1.0),
      73          36 :   cycleMaxPsfFraction_p(0.8),
      74          36 :   donePSF_p(false),
      75          36 :   modified_p(true),
      76          36 :   dataPolRep_p(StokesImageUtil::CIRCULAR),
      77          36 :   workDirOnNFS_p(false), useMem_p(false), tileVol_p(1000000)
      78             :  {
      79             :    
      80          36 :  }
      81             : 
      82           0 : void ImageSkyModel::setMaxNumberModels(const Int maxNumModels) {
      83           0 :   maxnmodels_p = maxNumModels;
      84           0 : }
      85             : 
      86           7 : Bool ImageSkyModel::add(ComponentList& compList)
      87             : {
      88           7 :   if(componentList_p==0) {
      89           7 :     componentList_p=new ComponentList(compList);
      90           7 :     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]=&image;
     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          32 : Int ImageSkyModel::add(ImageInterface<Float>& image, const Int maxNumXfr)
     119             : {
     120          32 :   Int thismodel=nmodels_p;
     121          32 :   nmodels_p++;
     122             :   /////Avoid using nfs file if memory is available
     123             :   try{
     124          32 :     Directory eldir(image.name());
     125          32 :     workDirOnNFS_p=eldir.isNFSMounted();
     126             :     //cerr << "Using NFS " << workDirOnNFS_p << endl;
     127          32 :   }
     128           0 :   catch(...){
     129           0 :     workDirOnNFS_p=false;
     130           0 :   }
     131             : 
     132          32 :   if(nmodels_p>maxnmodels_p) maxnmodels_p=nmodels_p;
     133             : 
     134          32 :   maxNumXFR_p=maxNumXfr;
     135             : 
     136          32 :   image_p.resize(nmodels_p); 
     137          32 :   cimage_p.resize(nmodels_p);
     138          32 :   cxfr_p.resize(nmodels_p*maxNumXFR_p);
     139          32 :   residual_p.resize(nmodels_p);
     140          32 :   residualImage_p.resize(nmodels_p);
     141          32 :   gS_p.resize(nmodels_p);
     142          32 :   psf_p.resize(nmodels_p);
     143          32 :   ggS_p.resize(nmodels_p);
     144          32 :   fluxScale_p.resize(nmodels_p);
     145          32 :   work_p.resize(nmodels_p);
     146          32 :   deltaimage_p.resize(nmodels_p);
     147          32 :   solve_p.resize(nmodels_p);
     148          32 :   doFluxScale_p.resize(nmodels_p);
     149          32 :   weight_p.resize(nmodels_p);
     150          32 :   beam_p.resize(nmodels_p);
     151             :   
     152          32 :   image_p[thismodel]=0;
     153          32 :   cimage_p[thismodel]=0;
     154          32 :   residual_p[thismodel]=0;
     155             : 
     156        1860 :   for (Int numXFR=0;numXFR<maxNumXFR_p;numXFR++) {
     157        1828 :     cxfr_p[thismodel*maxNumXFR_p+numXFR]=0;
     158             :   }
     159          32 :   residualImage_p[thismodel]=0;
     160          32 :   gS_p[thismodel]=0;
     161          32 :   psf_p[thismodel]=0;
     162          32 :   ggS_p[thismodel]=0;
     163          32 :   fluxScale_p[thismodel]=0;
     164          32 :   work_p[thismodel]=0;
     165          32 :   deltaimage_p[thismodel]=0;
     166          32 :   solve_p[thismodel]=true;
     167          32 :   doFluxScale_p[thismodel]=false;
     168          32 :   weight_p[thismodel]=0;
     169          32 :   beam_p[thismodel]=0;
     170             :   
     171             :   // Initialize image
     172          32 :   image_p[thismodel]=&image;
     173          32 :   AlwaysAssert(image_p[thismodel], AipsError);
     174          32 :   image_p[thismodel]->setUnits(Unit("Jy/pixel"));
     175          32 :   donePSF_p=false;
     176          32 :   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           0 : }
     192             :   
     193           0 :   ImageSkyModel::ImageSkyModel(const ImageSkyModel& other) : SkyModel() {
     194           0 :   operator=(other);
     195           0 : };
     196             : 
     197          36 : ImageSkyModel::~ImageSkyModel() {
     198          36 :   if(componentList_p) delete componentList_p; componentList_p=0;
     199          68 :   for (Int thismodel=0;thismodel<nmodels_p;thismodel++) {
     200          32 :     if(residualImage_p[thismodel]) {    
     201             :       //(residualImage_p[thismodel])->clearCache();
     202             :       //(residualImage_p[thismodel])->tempClose();
     203           0 :       delete residualImage_p[thismodel];       
     204             :     }
     205          32 :     residualImage_p[thismodel]=0;
     206          32 :     if(cimage_p[thismodel]){
     207             :       //(cimage_p[thismodel])->clearCache();
     208             :       //(cimage_p[thismodel])->tempClose();
     209          32 :       delete cimage_p[thismodel]; 
     210             :     }
     211          32 :     cimage_p[thismodel]=0;
     212        1860 :     for(Int numXFR=0;numXFR<maxNumXFR_p;numXFR++) {
     213        1828 :       if(cxfr_p[thismodel*maxNumXFR_p+numXFR])
     214           0 :         delete cxfr_p[thismodel*maxNumXFR_p+numXFR];
     215        1828 :       cxfr_p[thismodel*maxNumXFR_p+numXFR]=0;
     216             :     }
     217          32 :     if(gS_p[thismodel]) delete gS_p[thismodel]; gS_p[thismodel]=0;
     218          32 :     if(psf_p[thismodel]) delete psf_p[thismodel]; psf_p[thismodel]=0;
     219          32 :     if(ggS_p[thismodel]) delete ggS_p[thismodel]; ggS_p[thismodel]=0;
     220          32 :     if(fluxScale_p[thismodel]) delete fluxScale_p[thismodel]; fluxScale_p[thismodel]=0;
     221          32 :     if(work_p[thismodel]) delete work_p[thismodel]; work_p[thismodel]=0;
     222          32 :     if(deltaimage_p[thismodel]) delete deltaimage_p[thismodel]; deltaimage_p[thismodel]=0;
     223          32 :     if(weight_p[thismodel]) delete weight_p[thismodel]; weight_p[thismodel]=0;
     224          32 :     if(beam_p[thismodel]) delete beam_p[thismodel]; beam_p[thismodel]=0;
     225             :   }
     226          36 : };
     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           0 : }
     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          36 : Bool ImageSkyModel::isEmpty(Int model) 
     342             : {
     343          36 :   const IPosition tileShape = image(model).niceCursorShape();
     344          36 :   TiledLineStepper ls(image(model).shape(), tileShape, 0);
     345          36 :   RO_LatticeIterator<Float> li(image(model), ls);
     346             :   
     347          36 :   for(li.reset();!li.atEnd();li++) {
     348          36 :     if(max(abs(li.cursor()))>0.0) return false;
     349             :   }
     350           0 :   return true;
     351          36 : }
     352             : 
     353         170 : ImageInterface<Float>& ImageSkyModel::image(Int model) 
     354             : {
     355         170 :   AlwaysAssert(nmodels_p>0, AipsError);
     356         170 :   AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
     357         170 :   AlwaysAssert(image_p[model], AipsError);
     358         170 :   return *image_p[model];
     359             : };
     360             : 
     361             : 
     362             : 
     363         240 : ImageInterface<Complex>& ImageSkyModel::cImage(Int model) 
     364             : {
     365         240 :   AlwaysAssert(nmodels_p>0, AipsError);
     366         240 :   AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
     367             : 
     368             :   //if(model>0&&(cimage_p[model-1])) cimage_p[model-1]->tempClose();
     369             : 
     370         240 :   Double memoryMB=HostInfo::memoryFree()/1024/(MEMFACTOR*maxnmodels_p);
     371         240 :   if(cimage_p[model]==0) {
     372          32 :     Vector<Int> whichStokes(0);
     373          32 :     IPosition cimageShape;
     374          32 :     cimageShape=image_p[model]->shape();
     375             : 
     376             :     CoordinateSystem cimageCoord =
     377             :       StokesImageUtil::CStokesCoord(//cimageShape,
     378          64 :                                     image_p[model]->coordinates(),
     379             :                                     whichStokes,
     380          32 :                                     dataPolRep_p);
     381             :     //                              StokesImageUtil::CIRCULAR);
     382             : 
     383             :     /* STOKESDBG */ //cout << "ImageSkyModel::CImage : Correlation Planes 'whichStokes' : " << whichStokes << endl;
     384          32 :     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          64 :         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          32 :                               cimageCoord, (workDirOnNFS_p|| useMem_p) ? memoryMB : 0);
     397          32 :       cimagePtr->set(0.0);
     398          32 :       cimagePtr->setMaximumCacheSize(min(cimagePtr->shape().product()/10, 1000000));
     399             :     // cimagePtr->setMaximumCacheSize(cacheSize(model));
     400          32 :     AlwaysAssert(cimagePtr, AipsError);
     401          32 :     cimage_p[model] = cimagePtr;
     402          32 :     cimage_p[model]->setMiscInfo(image_p[model]->miscInfo());
     403          32 :     cimagePtr->clearCache();
     404          32 :   }
     405         240 :   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           0 :     } 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           0 : ImageInterface<Float>& ImageSkyModel::deltaImage(Int model) 
     633             : {
     634           0 :   AlwaysAssert(nmodels_p>0, AipsError);
     635           0 :   AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
     636             : 
     637             :   // if(model>0&&(deltaimage_p[model-1])) deltaimage_p[model-1]->tempClose();
     638             : 
     639           0 :   if(deltaimage_p[model]==0) {
     640           0 :     Double memoryMB=HostInfo::memoryTotal(true)/1024/(MEMFACTOR*maxnmodels_p);
     641             :         TempImage<Float>* deltaimagePtr = 
     642           0 :           new TempImage<Float> (TiledShape(image_p[model]->shape(),tileShape(model)), 
     643           0 :                             image_p[model]->coordinates(),
     644           0 :                                 (workDirOnNFS_p || useMem_p) ? memoryMB : 0);
     645             :     
     646           0 :     AlwaysAssert(deltaimagePtr, AipsError);
     647             :     //deltaimagePtr->setMaximumCacheSize(deltaimagePtr->shape().product()/10);
     648           0 :     deltaimagePtr->setMaximumCacheSize(cacheSize(model));
     649           0 :     deltaimagePtr->set(0.0);
     650           0 :     deltaimagePtr->clearCache();
     651           0 :     deltaimage_p[model] = deltaimagePtr;
     652             :   }
     653           0 :   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          36 : Bool ImageSkyModel::doFluxScale(Int model) 
     710             : {
     711          36 :   AlwaysAssert(nmodels_p>0, AipsError);
     712          36 :   AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
     713          36 :   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          22 : Bool ImageSkyModel::hasComponentList()
     724             : {
     725          22 :   return (componentList_p);
     726             : }
     727             : 
     728         236 : ComponentList& ImageSkyModel::componentList() 
     729             : {
     730         236 :   AlwaysAssert(componentList_p, AipsError);
     731         236 :   return *componentList_p;
     732             : }
     733             : 
     734          56 : Long ImageSkyModel::cacheSize(Int model){
     735          56 :   AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
     736             :   //  Long cachesize=(image_p[model]->shape()[0])*(image_p[model]->shape()[1])/4;
     737          56 :   Long cachesize=(image_p[model]->shape().product());
     738             :   // return Long((HostInfo::memoryFree()/(sizeof(Float)*16)*1024));
     739             : 
     740          56 :   return cachesize;
     741             : }
     742             : 
     743          88 : IPosition ImageSkyModel::tileShape(Int model){
     744          88 :   AlwaysAssert((model>-1)&&(model<nmodels_p), AipsError);
     745          88 :   Int x=Int(sqrt(Double(tileVol_p)));
     746         176 :   IPosition shp(4, min(image_p[model]->shape()(0), x), min(image_p[model]->shape()(1), x), 1, 1);
     747             :   
     748          88 :   return shp;
     749             : }
     750             : 
     751          14 : void ImageSkyModel::setMemoryUse(Bool useMem){
     752          14 :   useMem_p=useMem;
     753          14 : }
     754           0 : void ImageSkyModel::setTileVol(Int  tileVol){
     755           0 :   tileVol_p=tileVol;
     756           0 : }
     757             : 
     758             : } //# NAMESPACE CASA - END
     759             : 

Generated by: LCOV version 1.16