Line data    Source code

       1             : //# ClarkCleanAlgorithm.cc: implementation of ClarkCleanAlgorithm.h
       2             : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2003
       3             : //# Associated Universities, Inc. Washington DC, USA.
       4             : //#
       5             : //# This library is free software; you can redistribute it and/or modify it
       6             : //# under the terms of the GNU Library General Public License as published by
       7             : //# the Free Software Foundation; either version 2 of the License, or (at your
       8             : //# option) any later version.
       9             : //#
      10             : //# This library is distributed in the hope that it will be useful, but WITHOUT
      11             : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
      12             : //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
      13             : //# License for more details.
      14             : //#
      15             : //# You should have received a copy of the GNU Library General Public License
      16             : //# along with this library; if not, write to the Free Software Foundation,
      17             : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
      18             : //#
      19             : //# Correspondence concerning AIPS++ should be addressed as follows:
      20             : //#        Internet email: 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             : //# Includes
      29             : #include <casacore/casa/BasicSL/String.h>
      30             : #include <casacore/casa/Arrays/Array.h>
      31             : #include <casacore/lattices/Lattices/PagedArray.h>
      32             : #include <casacore/casa/Arrays/ArrayMath.h>
      33             : #include <casacore/casa/Arrays/Array.h>
      34             : #include <casacore/lattices/Lattices/ArrayLattice.h>
      35             : #include <casacore/lattices/Lattices/PagedArray.h>
      36             : #include <casacore/lattices/LEL/LatticeExpr.h>
      37             : #include <casacore/lattices/LEL/LatticeExprNode.h>
      38             : #include <casacore/casa/OS/File.h>
      39             : #include <casacore/casa/Exceptions/Error.h>
      40             : #include <casacore/casa/BasicSL/String.h>
      41             : #include <casacore/casa/Utilities/Assert.h>
      42             : #include <casacore/casa/OS/HostInfo.h>
      43             : 
      44             : #include <sstream>
      45             : #include <casacore/casa/Logging/LogMessage.h>
      46             : #include <casacore/casa/Logging/LogIO.h>
      47             : #include <casacore/casa/Logging/LogSink.h>
      48             : 
      49             : #include <synthesis/MeasurementEquations/LatConvEquation.h>
      50             : #include <synthesis/MeasurementEquations/ClarkCleanLatModel.h>
      51             : 
      52             : #include <synthesis/MeasurementComponents/ClarkCleanAlgorithm.h>
      53             : 
      54             : using namespace casacore;
      55             : namespace casa { //# NAMESPACE CASA - BEGIN
      56             : 
      57           0 : ClarkCleanAlgorithm::ClarkCleanAlgorithm() : model_sl_p(0), 
      58           0 :   myName("Clark Clean")
      59             : {
      60             : // Default constructor
      61             : //
      62           0 :   cache_p = HostInfo::memoryTotal(true)*1024/(16*16);
      63           0 : };
      64             : 
      65           0 : ClarkCleanAlgorithm::~ClarkCleanAlgorithm() 
      66             : {
      67             : // Default desctructor
      68             : //
      69           0 :   if (model_sl_p) delete model_sl_p;
      70           0 : };
      71             : 
      72           0 : void ClarkCleanAlgorithm::get() 
      73             : {
      74             : // Get the input data and parameters from the controller
      75             : //
      76             :   // Get the input parameters and data
      77           0 :   applicator.get(residual_sl);
      78           0 :   applicator.get(psf_sf);
      79           0 :   applicator.get(mask);
      80           0 :   applicator.get(gain);
      81           0 :   applicator.get(threshold);
      82           0 :   applicator.get(numberIterations);
      83           0 :   applicator.get(chan);
      84           0 :   applicator.get(nchan);
      85           0 :   if (model_sl_p) delete model_sl_p;
      86           0 :   model_sl_p = new PagedArray<Float>(TiledShape(residual_sl.shape()));
      87           0 :   model_sl_p->set(0.0f);
      88           0 :   model_sl_p->setMaximumCacheSize(cache_p);
      89           0 : }
      90             : 
      91           0 : void ClarkCleanAlgorithm::put() 
      92             : {
      93             : // Return the results to the controller
      94             : //
      95           0 :   applicator.put(model_sl_p->get());
      96           0 : }
      97             : 
      98           0 : String& ClarkCleanAlgorithm::name()
      99             : {
     100             : // Return the algorithm name
     101             : //
     102           0 : return myName;
     103             : };
     104             : 
     105           0 : void ClarkCleanAlgorithm::task()
     106             : {
     107             : // Do the parallelized part of the Clark CLEAN, acting on
     108             : // the local data obtained from the controller.
     109             : //
     110           0 :   LogIO os(LogOrigin("task","solve in parallel",WHERE));
     111           0 :   PagedArray<Float> dirty_sl(TiledShape(residual_sl.shape()));
     112           0 :   dirty_sl.setMaximumCacheSize(cache_p);
     113           0 :   dirty_sl.put(residual_sl);
     114           0 :   PagedArray<Float> resid_sl(TiledShape(residual_sl.shape()));
     115           0 :   resid_sl.setMaximumCacheSize(cache_p);
     116           0 :   resid_sl.put(residual_sl);
     117             : 
     118           0 :   ArrayLattice<Float> al_psf_sl(psf_sf);
     119             :   Float psfmax;
     120             :   {
     121           0 :     LatticeExprNode node = max(al_psf_sl);
     122           0 :     psfmax = node.getFloat();
     123             :   }
     124           0 :   if(nchan>1) {
     125           0 :     os<<"Processing channel "<<chan+1<<" of "<<nchan<<LogIO::POST;
     126             :   }
     127           0 :   if(psfmax==0.0) {
     128           0 :     os << "No data for this channel: skipping" << LogIO::POST;
     129             :   } else {
     130           0 :     LatConvEquation eqn(al_psf_sl, dirty_sl);
     131           0 :     ClarkCleanLatModel cleaner(*model_sl_p);
     132           0 :     ArrayLattice<Float> latMask(mask);
     133           0 :     if (mask.nelements() > 1) {
     134           0 :       cleaner.setMask(latMask);
     135             :     };
     136             :     
     137           0 :     cleaner.setGain(gain);
     138           0 :     cleaner.setNumberIterations(numberIterations);
     139           0 :     cleaner.setThreshold(threshold);
     140           0 :     cleaner.setPsfPatchSize(IPosition(2,51)); 
     141           0 :     cleaner.setHistLength(1024);
     142           0 :     cleaner.setMaxNumPix(32*1024);
     143           0 :     cleaner.solve(eqn);
     144           0 :     cleaner.setChoose(false);
     145             :     os << LogIO::NORMAL    // Loglevel INFO
     146             :        << "Clean used " << cleaner.numberIterations() << " iterations" 
     147           0 :        << " to get to a max residual of " << cleaner.threshold() 
     148           0 :        << LogIO::POST;
     149             :     // cleaner.getModel(image);
     150           0 :     eqn.residual(resid_sl, cleaner);
     151             :   }
     152           0 : };
     153             : 
     154             : 
     155             : 
     156             : 
     157             : } //# NAMESPACE CASA - END
     158             :